source: tags/ORCHIDEE/src_sechiba/sechiba.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: 68.2 KB
Line 
1!! This module computes continental processes SECHIBA
2!!
3!! See also this graph <a href="http:shema.gif"><img SRC="  http:shema.gif" WIDTH="25" LENGTH="25"></A>
4!!
5!! @author Marie-Alice Foujols and Jan Polcher
6!! @Version : $Revision: 1.46 $, $Date: 2010/05/07 08:28:13 $
7!!
8!! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/sechiba.f90,v 1.46 2010/05/07 08:28:13 ssipsl Exp $
9!! IPSL (2006)
10!!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
11!!
12MODULE sechiba
13
14  ! routines called : diffuco_main, enerbil_main, hydrolc_main, enrbil_fusion, condveg_main, thermosoil_main
15  !
16  USE ioipsl
17  !
18  ! modules used :
19  USE constantes
20  USE constantes_veg
21  USE constantes_co2
22  USE diffuco
23  USE condveg
24  USE enerbil
25  USE hydrol
26  USE hydrolc
27  USE thermosoil
28  USE sechiba_io
29  USE slowproc
30  USE routing
31!  USE write_field_p, only : WriteFieldI_p
32
33
34  IMPLICIT NONE
35
36  ! public routines :
37  ! sechiba_main
38  ! sechiba_clear
39  PRIVATE
40  PUBLIC sechiba_main,sechiba_clear
41
42  ! Index arrays we need internaly
43  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexveg       !! indexing array for the 3D fields of vegetation
44  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnobio     !! indexing array for the 3D fields of other surfaces (ice, lakes, ...)
45  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsoil      !! indexing array for the 3D fields of soil types
46  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexgrnd      !! indexing array for the 3D ground heat profiles
47  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlayer     !! indexing array for the 3D fields of soil layers
48  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexalb       !! indexing array for the 2 fields of albedo
49
50  ! three dimensions array allocated, computed, saved and got in sechiba module
51
52  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: assim_param    !! min+max+opt temps, vcmax, vjmax for photosynthesis
53
54  ! two dimensions array allocated, computed, saved and got in sechiba module
55
56  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: veget          !! Fraction of vegetation type       
57  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: veget_max      !! Max. fraction of vegetation type (LAI -> infty)
58  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: height         !! Vegetation Height (m)
59  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lai            !! Surface foliere
60  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: gpp            !! STOMATE: GPP. gC/m**2 of total area
61  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: humrel         !! Relative humidity
62  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vegstress      !! Vegetation moisture stress (only for vegetation growth)
63  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: soiltype       !! Map of soil types, created in slowproc in the
64  !! order : silt, sand and clay.
65
66  !
67  ! one dimension array allocated, computed and used in sechiba module and passed to other
68  ! modules called
69
70  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta1         !! Snow resistance
71  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta4         !! Bare soil resistance
72  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: soilcap        !! Soil calorific capacity
73  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: soilflx        !! Soil flux
74  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: temp_sol       !! Soil temperature
75  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: qsurf          !! near soil air moisture
76  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: snow           !! Snow mass [Kg/m^2]
77  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: snow_age       !! Snow age
78  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: drysoil_frac   !! Fraction of visibly (albedo) Dry soil (Between 0 and 1)
79  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: rsol           !! resistance to bare soil evaporation
80  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evap_bare_lim  !! Bare soil stress
81  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evapot         !! Soil Potential Evaporation
82  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evapot_corr    !! Soil Potential Evaporation Correction (Milly 1992)
83  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapsno       !! Snow evaporation
84  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapnu        !! Bare soil evaporation
85  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: t2mdiag        !! 2 meter temperature
86  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: tot_melt       !! Total melt
87  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta          !! Resistance coefficient
88  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: valpha         !! Resistance coefficient
89  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: fusion         !!
90  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: rau            !! Density
91  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: deadleaf_cover !! Fraction of soil covered by dead leaves
92  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: co2_flux       !! CO2 flux (gC/m**2 of average ground/s)
93  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: totfrac_nobio  !! Total fraction of continental ice+lakes+cities+...
94  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: runoff         !! Surface runoff generated by hydrol
95  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: drainage       !! Deep drainage generated by hydrol
96  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: returnflow     !! Routed water which returns into the soil
97  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: irrigation     !! irrigation going back into the soils
98  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: emis           !! Surface emissivity
99  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: z0             !! Surface roughness
100  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: roughheight    !! Effective height for roughness
101  !
102  ! Arrays which are diagnostics from the physical processes for
103  ! for the biological processes. They are not saved in the restart file because at the first time step,
104  ! they are recalculated. However, they must be saved in memory as they are in slowproc which is called
105  ! before the modules which calculate them.
106  !
107  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: shumdiag     !! Relative soil moisture
108  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: litterhumdiag!! litter humidity
109  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: stempdiag    !! Temperature which controls canopy evolution
110
111  ! two dimensions array allocated, computed and used in sechiba module and passed to other
112  ! modules called
113
114  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: qsintveg       !! Water on vegetation due to interception
115  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta2         !! Interception resistance
116  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta3         !! Vegetation resistance
117  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbetaco2       !! STOMATE: Vegetation resistance to CO2
118  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: cimean         !! STOMATE: Mean ci
119  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vevapwet       !! Interception
120  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: transpir       !! Transpiration
121  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: qsintmax       !! Maximum water on vegetation for interception
122  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: rveget         !! Vegetation resistance
123  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: snow_nobio     !! Water balance over other surface types (that is snow !)
124  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: snow_nobio_age !! Snow age on other surface types
125  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: frac_nobio     !! Fraction of continental ice, lakes, ...
126  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: albedo         !! Surface albedo for visible and NIR
127  !
128  ! variables used inside sechiba module : declaration and initialisation
129  !
130  LOGICAL, SAVE                                   :: l_first_sechiba = .TRUE.!! Initialisation has to be done one time
131  CHARACTER(LEN=80) , SAVE                             :: var_name                !! To store variables names for I/O
132
133  LOGICAL, SAVE                                   :: river_routing           !! Flag that decides if we route.
134  LOGICAL, SAVE                                   :: hydrol_cwrr             !! Selects the CWRR hydrology.
135
136  LOGICAL, SAVE                                   :: myfalse=.FALSE. 
137  LOGICAL, SAVE                                   :: mytrue=.TRUE.
138
139CONTAINS
140
141  !! Main routine for *sechiba* module.
142  !!
143  !! Should be called two times:
144  !! - first time to have initial values
145  !! - second time to have complete algorithm
146  !!
147  !! Algorithm:
148  !! 3 series of called SECHIBA processes
149  !! - initialisation (first time only)
150  !! - time step evolution (every time step)
151  !! - prepares output and storage of restart arrays (last time only)
152  !!
153  !! One serie consists of:
154  !! - call sechiba_var_init to do some initialisation
155  !! - call slowproc_main to do some daily initialisation
156  !! - call diffuco_main for diffusion coefficient calculation
157  !! - call enerbil_main for energy bilan calculation
158  !! - call hydrolc_main for hydrologic processes calculation
159  !! - call enerbil_fusion : last part with fusion
160  !! - call condveg_main for surface conditions such as roughness, albedo, and emmisivity
161  !! - call thermosoil_main for soil thermodynamic calculation
162  !! - call sechiba_end to swap new fields to previous
163  !!
164  !! @call sechiba_var_init
165  !! @call slowproc_main
166  !! @call diffuco_main
167  !! @call enerbil_main
168  !! @call hydrolc_main
169  !! @call enerbil_fusion
170  !! @call condveg_main
171  !! @call thermosoil_main
172  !! @call sechiba_end
173  !!
174  SUBROUTINE sechiba_main (kjit, kjpij, kjpindex, index, dtradia, date0, &
175       & ldrestart_read, ldrestart_write, control_in, &
176       & lalo, contfrac, neighbours, resolution,&
177       ! First level conditions
178! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
179!     & zlev, u, v, qair, temp_air, epot_air, ccanopy, &
180    & zlev, u, v, qair, q2m, t2m, temp_air, epot_air, ccanopy, &
181         ! Variables for the implicit coupling
182    & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
183         ! Rain, snow, radiation and surface pressure
184    & precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
185         ! Output : Fluxes
186    & vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
187         ! Surface temperatures and surface properties
188    & tsol_rad, temp_sol_new, qsurf_out, albedo_out, emis_out, z0_out, &
189         ! File ids
190    & rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
191
192    ! interface description for dummy arguments
193    ! input scalar
194    INTEGER(i_std), INTENT(in)                               :: kjit             !! Time step number
195    INTEGER(i_std), INTENT(in)                               :: kjpij            !! Total size of size. This is the un-compressed grid
196    INTEGER(i_std), INTENT(in)                               :: kjpindex         !! Domain size
197    INTEGER(i_std),INTENT (in)                               :: rest_id          !! _Restart_ file identifier
198    INTEGER(i_std),INTENT (in)                               :: hist_id          !! _History_ file identifier
199    INTEGER(i_std),INTENT (in)                               :: hist2_id         !! _History_ file 2 identifier
200    INTEGER(i_std),INTENT (in)                               :: rest_id_stom     !! STOMATE's _Restart_ file identifier
201    INTEGER(i_std),INTENT (in)                               :: hist_id_stom     !! STOMATE's _History_ file identifier
202    INTEGER(i_std),INTENT(in)                                :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file identifier
203    REAL(r_std), INTENT (in)                                 :: dtradia          !! Time step in seconds
204    REAL(r_std), INTENT (in)                                 :: date0            !! initial date
205    LOGICAL, INTENT(in)                                     :: ldrestart_read   !! Logical for _restart_ file to read
206    LOGICAL, INTENT(in)                                     :: ldrestart_write  !! Logical for _restart_ file to write
207    TYPE(control_type), INTENT(in)                          :: control_in       !! Flags that (de)activate parts of the model
208    ! input fields
209    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)          :: lalo             !! Geographical coordinates
210    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: contfrac         !! Fraction of continent in the grid
211    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index            !! Indeces of the points on the map
212    !
213    INTEGER(i_std), DIMENSION (kjpindex,8), INTENT(in)       :: neighbours       !! neighoring grid points if land
214    REAL(r_std), DIMENSION (kjpindex,2), INTENT(in)          :: resolution       !! size in x an y of the grid (m)
215    !
216    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u                !! Lowest level wind speed
217    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: v                !! Lowest level wind speed
218    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: zlev             !! Height of first layer
219    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair             !! Lowest level specific humidity
220! Ajout Nathalie - Juin 2006 - Q2M/t2m pour calcul Rveget
221    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q2m              !! 2m specific humidity
222    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: t2m              !! 2m air temperature
223    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: precip_rain      !! Rain precipitation
224    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: precip_snow      !! Snow precipitation
225    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: lwdown           !! Down-welling long-wave flux
226    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swnet            !! Net surface short-wave flux
227    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swdown           !! Down-welling surface short-wave flux
228    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_air         !! Air temperature in Kelvin
229    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: epot_air         !! Air potential energy
230    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ccanopy          !! CO2 concentration in the canopy
231    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petAcoef         !! Coeficients A from the PBL resolution
232    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqAcoef         !! One for T and another for q
233    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petBcoef         !! Coeficients B from the PBL resolution
234    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqBcoef         !! One for T and another for q
235    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: tq_cdrag         !! This is the cdrag without the wind multiplied
236    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: pb               !! Lowest level pressure
237    ! output fields
238    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: coastalflow      !! Diffuse water flow to the oceans
239    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: riverflow        !! River outflow to the oceans
240    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: tsol_rad         !! Radiative surface temperature
241    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapp           !! Total of evaporation
242    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: temp_sol_new     !! New soil temperature
243    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qsurf_out        !! Surface specific humidity
244    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0_out           !! Surface roughness (output diagnostic)
245    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out)         :: albedo_out       !! Albedo (output diagnostic)
246    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxsens         !! Sensible chaleur flux
247    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxlat          !! Latent chaleur flux
248    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: emis_out         !! Emissivity
249
250    REAL(r_std), ALLOCATABLE, DIMENSION (:)                  :: runoff1,drainage1, soilcap1,soilflx1
251    REAL(r_std), ALLOCATABLE, DIMENSION (:,:)                :: shumdiag1
252
253    REAL(r_std), DIMENSION(kjpindex)                   :: histvar          !! computations for history files
254
255    IF (long_print) WRITE(numout,*) ' kjpindex =',kjpindex
256
257    ! do SECHIBA'S first call initialisation
258
259    IF (l_first_sechiba) THEN
260
261       CALL sechiba_init (kjit, ldrestart_read, kjpij, kjpindex, index, rest_id, control_in, lalo)
262
263       ALLOCATE(runoff1 (kjpindex),drainage1 (kjpindex), soilcap1 (kjpindex),soilflx1 (kjpindex))
264       ALLOCATE(shumdiag1(kjpindex,nbdl))
265       !
266       ! computes initialisation of energy bilan
267       !
268       IF (ldrestart_read) THEN
269         
270          IF (long_print) WRITE (numout,*) ' we have to read a restart file for SECHIBA variables'
271          var_name='soilcap' ;
272          CALL ioconf_setatt('UNITS', '-')
273          CALL ioconf_setatt('LONG_NAME','Soil calorific capacity')
274          soilcap1=val_exp
275          IF ( ok_var(var_name) ) THEN
276             CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., soilcap1, "gather", nbp_glo, index_g)
277             IF (MINVAL(soilcap1) < MAXVAL(soilcap1) .OR. MAXVAL(soilcap1) < val_exp) THEN
278                soilcap(:) = soilcap1(:)
279             ENDIF
280          ENDIF
281          !
282          var_name='soilflx' ;
283          CALL ioconf_setatt('UNITS', '-')
284          CALL ioconf_setatt('LONG_NAME','Soil flux')
285          soilflx1=val_exp
286          IF ( ok_var(var_name) ) THEN
287             CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., soilflx1, "gather", nbp_glo, index_g)
288             IF (MINVAL(soilflx1) < MAXVAL(soilflx1)  .OR. MAXVAL(soilflx1) < val_exp) THEN
289                soilflx(:) = soilflx1(:)
290             ENDIF
291          ENDIF
292          !
293          var_name='shumdiag' ;
294          CALL ioconf_setatt('UNITS', '-')
295          CALL ioconf_setatt('LONG_NAME','Relative soil moisture')
296          shumdiag1=val_exp
297          IF ( ok_var(var_name) ) THEN
298             CALL restget_p (rest_id, var_name, nbp_glo, nbdl, 1, kjit, .TRUE., shumdiag1, "gather", nbp_glo, index_g)
299             IF (MINVAL(shumdiag1) < MAXVAL(shumdiag1) .OR. MAXVAL(shumdiag1) < val_exp) THEN
300                shumdiag(:,:) = shumdiag1(:,:)
301             ENDIF
302          ENDIF
303       ENDIF
304
305       !
306       ! computes slow variables
307       !
308       CALL slowproc_main (kjit, kjpij, kjpindex, dtradia, date0, &
309            ldrestart_read, ldrestart_write, control%ok_co2, control%ok_stomate, &
310            index, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
311            t2mdiag, t2mdiag, temp_sol, stempdiag, &
312            vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
313            deadleaf_cover, &
314            assim_param, &
315            lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
316            rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
317            co2_flux)
318       !
319       ! computes initialisation of diffusion coeff
320       !
321
322       CALL diffuco_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, u, v, &
323! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
324!           & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, pb ,  &
325            & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, q2m, t2m, pb ,  &
326            & rsol, evap_bare_lim, evapot, snow, frac_nobio, snow_nobio, totfrac_nobio, &
327            & swnet, swdown, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
328            & vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, rveget, cimean, rest_id, hist_id, hist2_id)
329       !
330       ! computes initialisation of energy bilan
331       !
332       CALL enerbil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, &
333            & index, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef,&
334            & qair, peqAcoef, peqBcoef, pb, rau, vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, &
335            & cimean, ccanopy, emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
336            & vevapp, transpir, gpp, vevapnu, vevapwet, vevapsno, t2mdiag, temp_sol, tsol_rad, &
337            & temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id) 
338       !
339       ! computes initialisation of hydrologie
340       !
341       IF (ldrestart_read) THEN
342         
343          var_name='runoff' ;
344          CALL ioconf_setatt('UNITS', 'mm/d')
345          CALL ioconf_setatt('LONG_NAME','Complete runoff')
346          runoff1=val_exp
347          IF ( ok_var(var_name) ) THEN
348             CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff1, "gather", nbp_glo, index_g)
349             IF (MINVAL(runoff1) < MAXVAL(runoff1) .OR. MAXVAL(runoff1) < val_exp) THEN
350                runoff(:) = runoff1(:)
351             ENDIF
352          ENDIF
353
354          var_name='drainage' ;
355          CALL ioconf_setatt('UNITS', 'mm/d')
356          CALL ioconf_setatt('LONG_NAME','Deep drainage')
357          drainage1=val_exp
358          IF ( ok_var(var_name) ) THEN
359             CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage1, "gather", nbp_glo, index_g)
360             IF (MINVAL(drainage1) < MAXVAL(drainage1) .OR. MAXVAL(drainage1) < val_exp) THEN
361                drainage(:) = drainage1(:)
362             ENDIF
363          ENDIF
364
365          IF ( ok_var("shumdiag") ) THEN
366             IF (MINVAL(shumdiag1) < MAXVAL(shumdiag1) .OR. MAXVAL(shumdiag1) < val_exp) THEN
367                shumdiag(:,:) = shumdiag1(:,:)
368             ENDIF
369          ENDIF
370       ENDIF
371       !
372       IF ( .NOT. hydrol_cwrr ) THEN
373          CALL hydrolc_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, &
374               & temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
375               & qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
376               & tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
377               & vegstress, rsol, drysoil_frac, evapot, evapot_corr, shumdiag, litterhumdiag, soilcap, rest_id, hist_id, hist2_id) 
378          evap_bare_lim(:) = -un
379       ELSE
380          CALL hydrol_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, indexsoil, indexlayer, &
381               & temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
382               & qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
383               & tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
384               & vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, &
385               & shumdiag, litterhumdiag, soilcap, soiltype, rest_id, hist_id, hist2_id)
386       ENDIF
387       IF (ldrestart_read) THEN
388         
389          IF ( ok_var("runoff") ) THEN
390             IF (MINVAL(runoff1) < MAXVAL(runoff1) .OR. MAXVAL(runoff1) < val_exp) THEN
391                runoff(:) = runoff1(:)
392             ENDIF
393          ENDIF
394
395          IF ( ok_var("drainage") ) THEN
396             IF (MINVAL(drainage1) < MAXVAL(drainage1) .OR. MAXVAL(drainage1) < val_exp) THEN
397                drainage(:) = drainage1(:)
398             ENDIF
399          ENDIF
400         
401          IF ( ok_var("shumdiag") ) THEN
402             IF (MINVAL(shumdiag1) < MAXVAL(shumdiag1) .OR. MAXVAL(shumdiag1) < val_exp) THEN
403                shumdiag(:,:) = shumdiag1(:,:)
404             ENDIF
405          ENDIF
406       ENDIF
407
408       !
409       ! computes initialisation of condveg
410       !
411       CALL condveg_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, &
412            & lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
413            & zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
414            & drysoil_frac, height, deadleaf_cover, emis, albedo, z0, roughheight, soiltype, rest_id, hist_id, hist2_id)
415       !
416       ! computes initialisation of Soil Thermodynamic
417       !
418       CALL thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, &
419            & index, indexgrnd, temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id) 
420
421       IF (ldrestart_read) THEN
422          IF ( ok_var("soilcap") ) THEN
423             IF (MINVAL(soilcap1) < MAXVAL(soilcap1) .OR. MAXVAL(soilcap1) < val_exp) THEN
424                soilcap(:) = soilcap1(:)
425             ENDIF
426          ENDIF
427          !
428          IF ( ok_var("soilflx") ) THEN
429             IF (MINVAL(soilflx1) < MAXVAL(soilflx1)  .OR. MAXVAL(soilflx1) < val_exp) THEN
430                soilflx(:) = soilflx1(:)
431             ENDIF
432          ENDIF
433       ENDIF
434
435       !
436       ! If we chose to route the water then we call the module. Else variables
437       ! are set to zero.
438       !
439       !
440       IF ( river_routing .AND. nbp_glo .GT. 1) THEN
441          CALL routing_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, &
442               & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, runoff, &
443               & drainage, evapot_corr, precip_rain, humrel, &
444               & stempdiag, returnflow, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
445       ELSE
446          riverflow(:) = zero
447          coastalflow(:) = zero
448          returnflow(:) = zero
449          irrigation(:) = zero
450       ENDIF
451       !
452       ! Write the internal variables into the output fields
453       !
454       z0_out(:) = z0(:)
455       emis_out(:) = emis(:)
456       albedo_out(:,:) = albedo(:,:) 
457       qsurf_out(:) = qsurf(:)
458
459       DEALLOCATE(runoff1,drainage1,soilcap1,soilflx1)
460       DEALLOCATE(shumdiag1)
461
462       !
463       ! This line should remain last as it ends the initialisation and returns to the caling
464       ! routine.
465       !
466       RETURN
467       !
468    ENDIF
469
470    !
471    ! computes some initialisation every SECHIBA's call
472    !
473    CALL sechiba_var_init (kjpindex, rau, pb, temp_air) 
474
475    !
476    ! computes diffusion coeff
477    !
478    CALL diffuco_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index, indexveg, u, v, &
479! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
480!        & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, pb ,  &
481         & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, q2m, t2m, pb ,  &
482         & rsol, evap_bare_lim, evapot, snow, frac_nobio, snow_nobio, totfrac_nobio, &
483         & swnet, swdown, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
484         & vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, rveget, cimean, rest_id, hist_id, hist2_id)
485
486    !
487    ! computes energy bilan
488    !
489
490    CALL enerbil_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, &
491         & index, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, &
492         & qair, peqAcoef, peqBcoef, pb, rau, vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, &
493         & cimean, ccanopy, emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
494         & vevapp, transpir, gpp, vevapnu, vevapwet, vevapsno, t2mdiag, temp_sol, tsol_rad, &
495         & temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id)
496    !
497    ! computes hydrologie
498    !
499    IF ( .NOT. hydrol_cwrr ) THEN
500       CALL hydrolc_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index, indexveg, &
501            & temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
502            & qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
503            & tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
504            & vegstress, rsol, drysoil_frac, evapot, evapot_corr, shumdiag, litterhumdiag, soilcap, rest_id, hist_id, hist2_id) 
505       evap_bare_lim(:) = -un
506    ELSE
507       CALL hydrol_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index, indexveg, indexsoil, indexlayer, &
508            & temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
509            & qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
510            & tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
511            & vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, &
512            & shumdiag, litterhumdiag, soilcap, soiltype, rest_id, hist_id, hist2_id)
513    ENDIF
514    !
515    ! computes last part of energy bilan
516    !
517    CALL enerbil_fusion (kjpindex, dtradia, tot_melt, soilcap, temp_sol_new, fusion)
518
519    !
520    ! computes condveg
521    !
522    CALL condveg_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index,&
523         & lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
524         & zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
525         & drysoil_frac, height, deadleaf_cover, emis, albedo, z0, roughheight, soiltype, rest_id, hist_id, hist2_id)
526    !
527    ! computes Soil Thermodynamic
528    !
529    CALL thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index, indexgrnd, &
530         & temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id) 
531
532    !
533    ! If we chose to route the water then we call the module. Else variables
534    ! are set to zero.
535    !
536    !
537    IF ( river_routing .AND. nbp_glo .GT. 1) THEN
538       CALL routing_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index, &
539            & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, runoff, &
540            & drainage, evapot_corr, precip_rain, humrel, &
541            & stempdiag, returnflow, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
542       !       returnflow(:) = returnflow(:) * 100.
543    ELSE
544       riverflow(:) = zero
545       coastalflow(:) = zero
546       returnflow(:) = zero
547       irrigation(:) = zero
548    ENDIF
549
550
551    !
552    ! computes slow variables
553    ! ok_co2 and ok_stomate are interpreted as flags that determine whether the
554    !   forcing files are written.
555    !
556    CALL slowproc_main (kjit, kjpij, kjpindex, dtradia, date0, &
557         ldrestart_read, myfalse, control%ok_co2, control%ok_stomate, &
558         index, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
559         t2mdiag, t2mdiag, temp_sol, stempdiag, &
560         vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
561         deadleaf_cover, &
562         assim_param, &
563         lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
564         rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
565         co2_flux)
566
567    !
568    ! call swap from new computed variables 
569    !
570    CALL sechiba_end (kjpindex, dtradia, temp_sol, temp_sol_new)
571    !
572    ! Write the internal variables into the output fields
573    !
574    z0_out(:) = z0(:)
575    emis_out(:) = emis(:)
576    albedo_out(:,:) = albedo(:,:) 
577    qsurf_out(:) = qsurf(:)
578    !
579    ! Writing the global variables on the history tape
580    !
581    !
582    IF ( .NOT. almaoutput ) THEN
583       CALL histwrite(hist_id, 'beta', kjit, vbeta, kjpindex, index)
584       CALL histwrite(hist_id, 'z0', kjit, z0, kjpindex, index)
585       CALL histwrite(hist_id, 'roughheight', kjit, roughheight, kjpindex, index)
586       CALL histwrite(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
587       CALL histwrite(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
588       CALL histwrite(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
589       CALL histwrite(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
590       CALL histwrite(hist_id, 'subli', kjit, vevapsno, kjpindex, index)
591       CALL histwrite(hist_id, 'evapnu', kjit, vevapnu, kjpindex, index)
592       CALL histwrite(hist_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
593       CALL histwrite(hist_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
594       CALL histwrite(hist_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
595       CALL histwrite(hist_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
596       CALL histwrite(hist_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
597       CALL histwrite(hist_id, 'vbeta4', kjit, vbeta4, kjpindex, index)   
598       CALL histwrite(hist_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
599       CALL histwrite(hist_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
600       CALL histwrite(hist_id, 'rsol', kjit, rsol, kjpindex, index)
601       CALL histwrite(hist_id, 'snow', kjit, snow, kjpindex, index)
602       CALL histwrite(hist_id, 'snowage', kjit, snow_age, kjpindex, index)
603       CALL histwrite(hist_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
604       CALL histwrite(hist_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
605       CALL histwrite(hist_id, 'soiltype',  kjit, soiltype, kjpindex*nstm, indexsoil)
606       IF ( control%ok_co2 ) THEN
607          CALL histwrite(hist_id, 'vbetaco2', kjit, vbetaco2, kjpindex*nvm, indexveg)   
608          CALL histwrite(hist_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
609          CALL histwrite(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)   
610       ENDIF
611       IF ( control%ok_stomate ) THEN
612          CALL histwrite(hist_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)   
613       ENDIF
614
615       histvar(:)=SUM(vevapwet(:,:),dim=2)/86400
616       CALL histwrite(hist_id, 'evspsblveg', kjit, histvar, kjpindex, index)
617
618       histvar(:)=(vevapnu(:)+vevapsno(:))/86400
619       CALL histwrite(hist_id, 'evspsblsoi', kjit, histvar, kjpindex, index)
620
621       histvar(:)=SUM(transpir(:,:),dim=2)/86400
622       CALL histwrite(hist_id, 'tran', kjit, histvar, kjpindex, index)
623
624       histvar(:)=SUM(veget_max(:,2:9),dim=2)*100*contfrac(:)
625       CALL histwrite(hist_id, 'treeFrac', kjit, histvar, kjpindex, index)
626
627       histvar(:)=SUM(veget_max(:,10:11),dim=2)*100*contfrac(:)
628       CALL histwrite(hist_id, 'grassFrac', kjit, histvar, kjpindex, index)
629
630       histvar(:)=SUM(veget_max(:,12:13),dim=2)*100*contfrac(:)
631       CALL histwrite(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)
632
633       histvar(:)=veget_max(:,1)*100*contfrac(:)
634       CALL histwrite(hist_id, 'baresoilFrac', kjit, histvar, kjpindex, index)
635
636       histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)*100*contfrac(:)
637       CALL histwrite(hist_id, 'residualFrac', kjit, histvar, kjpindex, index)
638
639    ELSE
640       CALL histwrite(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
641       CALL histwrite(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
642       CALL histwrite(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
643       CALL histwrite(hist_id, 'Qf', kjit, fusion, kjpindex, index)
644       CALL histwrite(hist_id, 'SWE', kjit, snow, kjpindex, index)
645       CALL histwrite(hist_id, 'ESoil', kjit, vevapnu, kjpindex, index)
646       CALL histwrite(hist_id, 'TVeg', kjit, transpir, kjpindex*nvm, indexveg)
647       CALL histwrite(hist_id, 'ECanop', kjit, vevapwet, kjpindex*nvm, indexveg)
648       CALL histwrite(hist_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
649       CALL histwrite(hist_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
650       IF ( control%ok_co2 ) THEN
651          CALL histwrite(hist_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
652       ENDIF
653       IF ( control%ok_stomate ) THEN
654             CALL histwrite(hist_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)   
655       ENDIF
656    ENDIF
657    !
658    IF ( hist2_id > 0 ) THEN
659       IF ( .NOT. almaoutput ) THEN
660          CALL histwrite(hist2_id, 'tsol_rad', kjit, tsol_rad, kjpindex, index)
661          CALL histwrite(hist2_id, 'qsurf', kjit, qsurf, kjpindex, index)
662          CALL histwrite(hist2_id, 'albedo', kjit, albedo, kjpindex*2, indexalb)
663          CALL histwrite(hist2_id, 'emis', kjit, emis, kjpindex, index)
664          !
665          CALL histwrite(hist2_id, 'beta', kjit, vbeta, kjpindex, index)
666          CALL histwrite(hist2_id, 'z0', kjit, z0, kjpindex, index)
667          CALL histwrite(hist2_id, 'roughheight', kjit, roughheight, kjpindex, index)
668          CALL histwrite(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
669          CALL histwrite(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
670          CALL histwrite(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
671          CALL histwrite(hist2_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
672          CALL histwrite(hist2_id, 'subli', kjit, vevapsno, kjpindex, index)
673          CALL histwrite(hist2_id, 'vevapnu', kjit, vevapnu, kjpindex, index)
674          CALL histwrite(hist2_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
675          CALL histwrite(hist2_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
676          CALL histwrite(hist2_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
677          CALL histwrite(hist2_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
678          CALL histwrite(hist2_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
679          CALL histwrite(hist2_id, 'vbeta4', kjit, vbeta4, kjpindex, index)   
680          CALL histwrite(hist2_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
681          CALL histwrite(hist2_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
682          CALL histwrite(hist2_id, 'rsol', kjit, rsol, kjpindex, index)
683          CALL histwrite(hist2_id, 'snow', kjit, snow, kjpindex, index)
684          CALL histwrite(hist2_id, 'snowage', kjit, snow_age, kjpindex, index)
685          CALL histwrite(hist2_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
686          CALL histwrite(hist2_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
687          CALL histwrite(hist2_id, 'soiltype',  kjit, soiltype, kjpindex*nstm, indexsoil)
688          IF ( control%ok_co2 ) THEN
689             CALL histwrite(hist2_id, 'vbetaco2', kjit, vbetaco2, kjpindex*nvm, indexveg)   
690             CALL histwrite(hist2_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
691             CALL histwrite(hist2_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)   
692          ENDIF
693          IF ( control%ok_stomate ) THEN
694             CALL histwrite(hist2_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)   
695          ENDIF
696       ELSE
697          CALL histwrite(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
698          CALL histwrite(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
699          CALL histwrite(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
700          CALL histwrite(hist2_id, 'Qf', kjit, fusion, kjpindex, index)
701          CALL histwrite(hist2_id, 'SWE', kjit, snow, kjpindex, index)
702          CALL histwrite(hist2_id, 'ESoil', kjit, vevapnu, kjpindex, index)
703          CALL histwrite(hist2_id, 'TVeg', kjit, transpir, kjpindex*nvm, indexveg)
704          CALL histwrite(hist2_id, 'ECanop', kjit, vevapwet, kjpindex*nvm, indexveg)
705          CALL histwrite(hist2_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
706          CALL histwrite(hist2_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
707          IF ( control%ok_co2 ) THEN
708             CALL histwrite(hist2_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
709          ENDIF
710          IF ( control%ok_stomate ) THEN
711             CALL histwrite(hist2_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)   
712          ENDIF
713       ENDIF
714    ENDIF
715    !
716    ! prepares restart file for the next simulation from SECHIBA and from other modules
717    !
718    IF (ldrestart_write) THEN
719
720       IF (long_print) WRITE (numout,*) ' we have to write a restart file '
721
722       !
723       ! call diffuco_main to write restart files
724       !
725
726       CALL diffuco_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, u, v, &
727! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
728!           & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, pb ,  &
729            & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, q2m, t2m, pb ,  &
730            & rsol, evap_bare_lim, evapot, snow, frac_nobio, snow_nobio, totfrac_nobio, &
731            & swnet, swdown, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
732            & vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, rveget, cimean, rest_id, hist_id, hist2_id)
733
734       !
735       ! call energy bilan to write restart files
736       !
737       CALL enerbil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, &
738            & index, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef,&
739            & qair, peqAcoef, peqBcoef, pb, rau, vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, &
740            & cimean, ccanopy, emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
741            & vevapp, transpir, gpp, vevapnu, vevapwet, vevapsno, t2mdiag, temp_sol, tsol_rad, &
742            & temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id) 
743
744       !
745       ! call hydrologie to write restart files
746       !
747       IF ( .NOT. hydrol_cwrr ) THEN
748          CALL hydrolc_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, &
749               & temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
750               & qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
751               & tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, &
752               & humrel, vegstress, rsol, drysoil_frac, evapot, evapot_corr, shumdiag, litterhumdiag, soilcap, &
753               & rest_id, hist_id, hist2_id)
754          evap_bare_lim(:) = -un
755       ELSE
756          CALL hydrol_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, indexsoil, indexlayer, &
757               & temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
758               & qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
759               & tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
760               & vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, &
761               & shumdiag, litterhumdiag, soilcap, soiltype, rest_id, hist_id, hist2_id)
762          rsol(:) = -un
763       ENDIF
764       !
765       ! call condveg to write restart files
766       !
767       CALL condveg_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, &
768            & lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
769            & zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
770            & drysoil_frac, height, deadleaf_cover, emis, albedo, z0, roughheight, soiltype, rest_id, hist_id, hist2_id)
771
772       !
773       ! call Soil Thermodynamic to write restart files
774       !
775       CALL thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexgrnd, &
776            & temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id) 
777
778       !
779       ! If we chose to route the water then we call the module. Else variables
780       ! are set to zero.
781       !
782       !
783       IF ( river_routing .AND. nbp_glo .GT. 1) THEN
784          CALL routing_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, &
785               & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, runoff, &
786               & drainage, evapot_corr, precip_rain, humrel, &
787               & stempdiag, returnflow, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
788       ELSE
789          riverflow(:) = zero
790          coastalflow(:) = zero
791          returnflow(:) = zero
792          irrigation(:) = zero
793       ENDIF
794
795       !
796       ! call slowproc_main to write restart files
797       !
798
799       CALL slowproc_main (kjit, kjpij, kjpindex, dtradia, date0, &
800            ldrestart_read, ldrestart_write, control%ok_co2, control%ok_stomate, &
801            index, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
802            t2mdiag, t2mdiag, temp_sol, stempdiag, &
803            vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
804            deadleaf_cover, &
805            assim_param, &
806            lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
807            rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
808            co2_flux)
809
810
811       var_name= 'shumdiag' 
812       CALL restput_p(rest_id, var_name, nbp_glo,   nbdl, 1, kjit,  shumdiag, 'scatter',  nbp_glo, index_g)
813
814       var_name= 'runoff' 
815       CALL restput_p(rest_id, var_name, nbp_glo,   1, 1, kjit,  runoff, 'scatter',  nbp_glo, index_g)
816       var_name= 'drainage' 
817       CALL restput_p(rest_id, var_name, nbp_glo,   1, 1, kjit,  drainage, 'scatter',  nbp_glo, index_g)
818    END IF
819
820    IF (long_print) WRITE (numout,*) ' sechiba_main done '
821
822  END SUBROUTINE sechiba_main
823
824  !! Initialisation for SECHIBA processes
825  !! - does dynamic allocation for local arrays
826  !! - reads _restart_ file or set initial values to a raisonable value
827  !! - reads initial map
828  !!
829  SUBROUTINE sechiba_init (kjit, ldrestart_read, kjpij, kjpindex, index, rest_id, control_in, lalo)
830
831    ! interface description
832    ! input scalar
833    INTEGER(i_std), INTENT (in)                         :: kjit               !! Time step number
834    LOGICAL,INTENT (in)                                :: ldrestart_read     !! Logical for restart file to read
835    INTEGER(i_std), INTENT (in)                         :: kjpij              !! Size of full domaine
836    INTEGER(i_std), INTENT (in)                         :: kjpindex           !! Domain size
837    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: index              !! Indeces of the points on the map
838    INTEGER(i_std), INTENT (in)                         :: rest_id            !! _Restart_ file identifier
839    TYPE(control_type), INTENT(in)                     :: control_in         !! Flags that (de)activate parts of the model
840    ! input fields
841    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)     :: lalo               !! Geographical coordinates
842    ! output scalar
843    ! output fields
844
845    ! local declaration
846    INTEGER(i_std)                                :: ier, ji, jv
847    !
848    ! initialisation
849    !
850    IF (l_first_sechiba) THEN
851       l_first_sechiba=.FALSE.
852    ELSE
853       WRITE (numout,*) ' l_first_sechiba false . we stop '
854       STOP 'sechiba_init'
855    ENDIF
856
857    ! 1. make dynamic allocation with good dimension
858
859    ! 1.0 The 3D vegetation indexation table
860
861    ALLOCATE (indexveg(kjpindex*nvm),stat=ier)
862    IF (ier.NE.0) THEN
863       WRITE (numout,*) ' error in indexveg allocation. We stop. We need kjpindex words = ',kjpindex*nvm
864       STOP 'sechiba_init'
865    END IF
866
867    ALLOCATE (indexsoil(kjpindex*nstm),stat=ier)
868    IF (ier.NE.0) THEN
869       WRITE (numout,*) ' error in indexsoil allocation. We stop. We need kjpindex words = ',kjpindex*nstm
870       STOP 'sechiba_init'
871    END IF
872
873    ALLOCATE (indexnobio(kjpindex*nnobio),stat=ier)
874    IF (ier.NE.0) THEN
875       WRITE (numout,*) ' error in indexnobio allocation. We stop. We need kjpindex words = ',kjpindex*nnobio
876       STOP 'sechiba_init'
877    END IF
878
879    ALLOCATE (indexgrnd(kjpindex*ngrnd),stat=ier)
880    IF (ier.NE.0) THEN
881       WRITE (numout,*) ' error in indexgrnd allocation. We stop. We need kjpindex words = ',kjpindex*ngrnd
882       STOP 'sechiba_init'
883    END IF
884
885    ALLOCATE (indexlayer(kjpindex*nslm),stat=ier)
886    IF (ier.NE.0) THEN
887       WRITE (numout,*) ' error in indexlayer allocation. We stop. We need kjpindex words = ',kjpindex*nslm
888       STOP 'sechiba_init'
889    END IF
890
891    ALLOCATE (indexalb(kjpindex*2),stat=ier)
892    IF (ier.NE.0) THEN
893       WRITE (numout,*) ' error in indexalb allocation. We stop. We need kjpindex words = ',kjpindex*2
894       STOP 'sechiba_init'
895    END IF
896
897    ! 1.1 one dimension array allocation with restartable value
898
899    IF (long_print) WRITE (numout,*) 'Allocation of 1D variables. We need for each kjpindex words = ',kjpindex
900    ALLOCATE (snow(kjpindex),stat=ier)
901    IF (ier.NE.0) THEN
902       WRITE (numout,*) ' error in snow allocation. We stop. We need kjpindex words = ',kjpindex
903       STOP 'sechiba_init'
904    END IF
905    snow(:) = undef_sechiba
906
907    ALLOCATE (snow_age(kjpindex),stat=ier)
908    IF (ier.NE.0) THEN
909       WRITE (numout,*) ' error in snow_age allocation. We stop. We need kjpindex words = ',kjpindex
910       STOP 'sechiba_init'
911    END IF
912    snow_age(:) = undef_sechiba
913
914    ALLOCATE (drysoil_frac(kjpindex),stat=ier)
915    IF (ier.NE.0) THEN
916       WRITE (numout,*) ' error in drysoil_frac allocation. We stop. We need kjpindex words = ',kjpindex
917       STOP 'sechiba_init'
918    END IF
919    drysoil_frac(:) = zero
920
921    ALLOCATE (rsol(kjpindex),stat=ier)
922    IF (ier.NE.0) THEN
923       WRITE (numout,*) ' error in rsol allocation. We stop. We need kjpindex words = ',kjpindex
924       STOP 'sechiba_init'
925    END IF
926
927    ALLOCATE (evap_bare_lim(kjpindex),stat=ier)
928    IF (ier.NE.0) THEN
929       WRITE (numout,*) ' error in evap_bare_lim allocation. We stop. We need kjpindex words = ',kjpindex
930       STOP 'sechiba_init'
931    END IF
932
933    ALLOCATE (evapot(kjpindex),stat=ier)
934    IF (ier.NE.0) THEN
935       WRITE (numout,*) ' error in evapot allocation. We stop. We need kjpindex words = ',kjpindex
936       STOP 'sechiba_init'
937    END IF
938    evapot(:) = undef_sechiba
939
940    ALLOCATE (evapot_corr(kjpindex),stat=ier)
941    IF (ier.NE.0) THEN
942       WRITE (numout,*) ' error in evapot_corr allocation. We stop. We need kjpindex words = ',kjpindex
943       STOP 'sechiba_init'
944    END IF
945
946    ALLOCATE (humrel(kjpindex,nvm),stat=ier)
947    IF (ier.NE.0) THEN
948       WRITE (numout,*) ' error in humrel allocation. We stop. we need kjpindex words = ',kjpindex
949       STOP 'sechiba_init'
950    END IF
951    humrel(:,:) = undef_sechiba
952
953    ALLOCATE (vegstress(kjpindex,nvm),stat=ier)
954    IF (ier.NE.0) THEN
955       WRITE (numout,*) ' error in vegstress allocation. We stop. we need kjpindex words = ',kjpindex
956       STOP 'sechiba_init'
957    END IF
958    vegstress(:,:) = undef_sechiba
959
960    ALLOCATE (soiltype(kjpindex,nstm),stat=ier)
961    IF (ier.NE.0) THEN
962       WRITE (numout,*) ' error in soiltype allocation. We stop. we need kjpindex words = ',kjpindex
963       STOP 'sechiba_init'
964    END IF
965    soiltype(:,:)=undef_sechiba
966
967    ALLOCATE (vbeta1(kjpindex),stat=ier)
968    IF (ier.NE.0) THEN
969       WRITE (numout,*) ' error in vbeta1 allocation. We stop. We need kjpindex words = ',kjpindex
970       STOP 'sechiba_init'
971    END IF
972
973    ALLOCATE (vbeta4(kjpindex),stat=ier)
974    IF (ier.NE.0) THEN
975       WRITE (numout,*) ' error in vbeta4 allocation. We stop. We need kjpindex words = ',kjpindex
976       STOP 'sechiba_init'
977    END IF
978
979    ALLOCATE (soilcap(kjpindex),stat=ier)
980    IF (ier.NE.0) THEN
981       WRITE (numout,*) ' error in soilcap allocation. We stop. We need kjpindex words = ',kjpindex
982       STOP 'sechiba_init'
983    END IF
984
985    ALLOCATE (soilflx(kjpindex),stat=ier)
986    IF (ier.NE.0) THEN
987       WRITE (numout,*) ' error in soilflx allocation. We stop. We need kjpindex words = ',kjpindex
988       STOP 'sechiba_init'
989    END IF
990
991    ALLOCATE (temp_sol(kjpindex),stat=ier)
992    IF (ier.NE.0) THEN
993       WRITE (numout,*) ' error in temp_sol allocation. We stop. We need kjpindex words = ',kjpindex
994       STOP 'sechiba_init'
995    END IF
996    temp_sol(:) = undef_sechiba
997
998    ALLOCATE (qsurf(kjpindex),stat=ier)
999    IF (ier.NE.0) THEN
1000       WRITE (numout,*) ' error in qsurf allocation. We stop. We need kjpindex words = ',kjpindex
1001       STOP 'sechiba_init'
1002    END IF
1003    qsurf(:) = undef_sechiba
1004
1005    ! 1.2 two dimensions array allocation with restartable value
1006
1007    ALLOCATE (qsintveg(kjpindex,nvm),stat=ier)
1008    IF (ier.NE.0) THEN
1009       WRITE (numout,*) ' error in qsintveg allocation. We stop. We need kjpindex x nvm words = ',&
1010            & kjpindex,' x ' ,nvm,' = ',kjpindex*nvm 
1011       STOP 'sechiba_init'
1012    END IF
1013    qsintveg(:,:) = undef_sechiba
1014
1015    ALLOCATE (vbeta2(kjpindex,nvm),stat=ier)
1016    IF (ier.NE.0) THEN
1017       WRITE (numout,*) ' error in vbeta2 allocation. We stop. We need kjpindex x nvm words = ',&
1018            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1019       STOP 'sechiba_init'
1020    END IF
1021
1022    ALLOCATE (vbeta3(kjpindex,nvm),stat=ier)
1023    IF (ier.NE.0) THEN
1024       WRITE (numout,*) ' error in vbeta3 allocation. We stop.We need kjpindex x nvm words = ',&
1025            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1026       STOP 'sechiba_init'
1027    END IF
1028
1029    ALLOCATE (vbetaco2(kjpindex,nvm),stat=ier)
1030    IF (ier.NE.0) THEN
1031       WRITE (numout,*) ' error in vbetaco2 allocation. We stop.We need kjpindex x nvm words = ',&
1032            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1033       STOP 'sechiba_init'
1034    END IF
1035
1036    ALLOCATE (cimean(kjpindex,nvm),stat=ier)
1037    IF (ier.NE.0) THEN
1038       WRITE (numout,*) ' error in cimean allocation. We stop.We need kjpindex x nvm words = ',&
1039            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1040       STOP 'sechiba_init'
1041    END IF
1042
1043    ALLOCATE (gpp(kjpindex,nvm),stat=ier)
1044    IF (ier.NE.0) THEN
1045       WRITE (numout,*) ' error in gpp allocation. We stop.We need kjpindex x nvm words = ',&
1046            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1047       STOP 'sechiba_init'
1048    END IF
1049    gpp(:,:) = undef_sechiba
1050
1051    ALLOCATE (veget(kjpindex,nvm),stat=ier)
1052    IF (ier.NE.0) THEN
1053       WRITE (numout,*) ' error in veget allocation. We stop. We need kjpindex x nvm words = ',&
1054            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1055       STOP 'sechiba_init'
1056    END IF
1057    veget(:,:)=undef_sechiba
1058
1059    ALLOCATE (veget_max(kjpindex,nvm),stat=ier)
1060    IF (ier.NE.0) THEN
1061       WRITE (numout,*) ' error in veget_max allocation. We stop. We need kjpindex x nvm words = ',&
1062            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1063       STOP 'sechiba_init'
1064    END IF
1065    veget_max(:,:)=undef_sechiba
1066
1067    ALLOCATE (lai(kjpindex,nvm),stat=ier)
1068    IF (ier.NE.0) THEN
1069       WRITE (numout,*) ' error in lai allocation. We stop. We need kjpindex x nvm words = ',&
1070            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1071       STOP 'sechiba_init'
1072    END IF
1073    lai(:,:)=undef_sechiba
1074
1075    ALLOCATE (height(kjpindex,nvm),stat=ier)
1076    IF (ier.NE.0) THEN
1077       WRITE (numout,*) ' error in height allocation. We stop. We need kjpindex x nvm words = ',&
1078            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1079       STOP 'sechiba_init'
1080    END IF
1081    height(:,:)=undef_sechiba
1082
1083    ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
1084    IF (ier.NE.0) THEN
1085       WRITE (numout,*) ' error in frac_nobio allocation. We stop. We need kjpindex words = ',kjpindex
1086       STOP 'sechiba_init'
1087    END IF
1088    frac_nobio(:,:) = undef_sechiba
1089
1090    ALLOCATE (albedo(kjpindex,2),stat=ier)
1091    IF (ier.NE.0) THEN
1092       WRITE (numout,*) ' error in albedo allocation. We stop. We need kjpindex words = ',kjpindex
1093       STOP 'sechiba_init'
1094    END IF
1095
1096    ALLOCATE (snow_nobio(kjpindex,nnobio),stat=ier)
1097    IF (ier.NE.0) THEN
1098       WRITE (numout,*) ' error in snow_nobio allocation. We stop. We need kjpindex words = ',kjpindex
1099       STOP 'sechiba_init'
1100    END IF
1101    snow_nobio(:,:) = undef_sechiba
1102
1103    ALLOCATE (snow_nobio_age(kjpindex,nnobio),stat=ier)
1104    IF (ier.NE.0) THEN
1105       WRITE (numout,*) ' error in snow_nobio_age allocation. We stop. We need kjpindex words = ',kjpindex
1106       STOP 'sechiba_init'
1107    END IF
1108    snow_nobio_age(:,:) = undef_sechiba
1109
1110    ALLOCATE (assim_param(kjpindex,nvm,npco2),stat=ier)
1111    IF (ier.NE.0) THEN
1112       WRITE (numout,*) ' error in assim_param allocation. We stop. We need kjpindex x nvm x npco2 words = ',&
1113            & kjpindex,' x ' ,nvm,' x ',npco2, ' = ',kjpindex*nvm*npco2
1114       STOP 'sechiba_init'
1115    END IF
1116
1117    ! 1.3 one dimension array allocation
1118
1119    ALLOCATE (vevapsno(kjpindex),stat=ier)
1120    IF (ier.NE.0) THEN
1121       WRITE (numout,*) ' error in vevapsno allocation. We stop. We need kjpindex words = ',kjpindex
1122       STOP 'sechiba_init'
1123    END IF
1124
1125    ALLOCATE (vevapnu(kjpindex),stat=ier)
1126    IF (ier.NE.0) THEN
1127       WRITE (numout,*) ' error in vevapnu allocation. We stop. We need kjpindex words = ',kjpindex
1128       STOP 'sechiba_init'
1129    END IF
1130
1131    ALLOCATE (t2mdiag(kjpindex),stat=ier)
1132    IF (ier.NE.0) THEN
1133       WRITE (numout,*) ' error in t2mdiag allocation. We stop. We need kjpindex words = ',kjpindex
1134       STOP 'sechiba_init'
1135    END IF
1136
1137    ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
1138    IF (ier.NE.0) THEN
1139       WRITE (numout,*) ' error in totfrac_nobio allocation. We stop. We need kjpindex words = ',kjpindex
1140       STOP 'sechiba_init'
1141    END IF
1142
1143    ALLOCATE (runoff(kjpindex),stat=ier)
1144    IF (ier.NE.0) THEN
1145       WRITE (numout,*) ' error in runoff allocation. We stop. We need kjpindex words = ',kjpindex
1146       STOP 'sechiba_init'
1147    END IF
1148
1149    ALLOCATE (drainage(kjpindex),stat=ier)
1150    IF (ier.NE.0) THEN
1151       WRITE (numout,*) ' error in drainage allocation. We stop. We need kjpindex words = ',kjpindex
1152       STOP 'sechiba_init'
1153    END IF
1154
1155    ALLOCATE (returnflow(kjpindex),stat=ier)
1156    IF (ier.NE.0) THEN
1157       WRITE (numout,*) ' error in returnflow allocation. We stop. We need kjpindex words = ',kjpindex
1158       STOP 'sechiba_init'
1159    END IF
1160    returnflow(:) = zero
1161
1162    ALLOCATE (irrigation(kjpindex),stat=ier)
1163    IF (ier.NE.0) THEN
1164       WRITE (numout,*) ' error in irrigation allocation. We stop. We need kjpindex words = ',kjpindex
1165       STOP 'sechiba_init'
1166    END IF
1167    irrigation(:) = zero
1168
1169    ALLOCATE (z0(kjpindex),stat=ier)
1170    IF (ier.NE.0) THEN
1171       WRITE (numout,*) ' error in z0 allocation. We stop. We need kjpindex words = ',kjpindex
1172       STOP 'sechiba_init'
1173    END IF
1174
1175    ALLOCATE (roughheight(kjpindex),stat=ier)
1176    IF (ier.NE.0) THEN
1177       WRITE (numout,*) ' error in roughheight allocation. We stop. We need kjpindex words = ',kjpindex
1178       STOP 'sechiba_init'
1179    END IF
1180
1181    ALLOCATE (emis(kjpindex),stat=ier)
1182    IF (ier.NE.0) THEN
1183       WRITE (numout,*) ' error in emis allocation. We stop. We need kjpindex words = ',kjpindex
1184       STOP 'sechiba_init'
1185    END IF
1186
1187    ALLOCATE (tot_melt(kjpindex),stat=ier)
1188    IF (ier.NE.0) THEN
1189       WRITE (numout,*) ' error in tot_melt allocation. We stop. We need kjpindex words = ',kjpindex
1190       STOP 'sechiba_init'
1191    END IF
1192
1193    ALLOCATE (valpha(kjpindex),stat=ier)
1194    IF (ier.NE.0) THEN
1195       WRITE (numout,*) ' error in valpha allocation. We stop. We need kjpindex words = ',kjpindex
1196       STOP 'sechiba_init'
1197    END IF
1198
1199    ALLOCATE (vbeta(kjpindex),stat=ier)
1200    IF (ier.NE.0) THEN
1201       WRITE (numout,*) ' error in vbeta allocation. We stop. We need kjpindex words = ',kjpindex
1202       STOP 'sechiba_init'
1203    END IF
1204
1205    ALLOCATE (fusion(kjpindex),stat=ier)
1206    IF (ier.NE.0) THEN
1207       WRITE (numout,*) ' error in fusion allocation. We stop. We need kjpindex words = ',kjpindex
1208       STOP 'sechiba_init'
1209    END IF
1210
1211    ALLOCATE (rau(kjpindex),stat=ier)
1212    IF (ier.NE.0) THEN
1213       WRITE (numout,*) ' error in rau allocation. We stop. We need kjpindex words = ',kjpindex
1214       STOP 'sechiba_init'
1215    END IF
1216
1217    ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
1218    IF (ier.NE.0) THEN
1219       WRITE (numout,*) ' error in deadleaf_cover allocation. We stop. We need kjpindex words = ',kjpindex
1220       STOP 'sechiba_init'
1221    END IF
1222
1223    ALLOCATE (stempdiag(kjpindex, nbdl),stat=ier)
1224    IF (ier.NE.0) THEN
1225       WRITE (numout,*) ' error in stempdiag allocation. We stop. We need kjpindex*nbdl words = ',&
1226            & kjpindex*nbdl
1227       STOP 'sechiba_init'
1228    END IF
1229    ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
1230    IF (ier.NE.0) THEN
1231       WRITE (numout,*) ' error in co2_flux allocation. We stop. We need kjpindex*nvm words = ' ,kjpindex*nvm
1232       STOP 'sechiba_init'
1233    END IF
1234    co2_flux(:,:)=zero
1235
1236    ALLOCATE (shumdiag(kjpindex,nbdl),stat=ier)
1237    IF (ier.NE.0) THEN
1238       WRITE (numout,*) ' error in shumdiag allocation. We stop. We need kjpindex*nbdl words = ',&
1239            & kjpindex*nbdl
1240       STOP 'sechiba_init'
1241    END IF
1242
1243    ALLOCATE (litterhumdiag(kjpindex),stat=ier)
1244    IF (ier.NE.0) THEN
1245       WRITE (numout,*) ' error in litterhumdiag allocation. We stop. We need kjpindex words = ',kjpindex
1246       STOP 'sechiba_init'
1247    END IF
1248
1249    ! 1.4 two dimensions array allocation
1250
1251    ALLOCATE (vevapwet(kjpindex,nvm),stat=ier)
1252    IF (ier.NE.0) THEN
1253       WRITE (numout,*) ' error in vevapwet allocation. We stop. We need kjpindex x nvm words = ',&
1254            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1255       STOP 'sechiba_init'
1256    END IF
1257
1258    ALLOCATE (transpir(kjpindex,nvm),stat=ier)
1259    IF (ier.NE.0) THEN
1260       WRITE (numout,*) ' error in transpir allocation. We stop. We need kjpindex x nvm words = ',&
1261            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1262       STOP 'sechiba_init'
1263    END IF
1264
1265    ALLOCATE (qsintmax(kjpindex,nvm),stat=ier)
1266    IF (ier.NE.0) THEN
1267       WRITE (numout,*) ' error in qsintmax allocation. We stop. We need kjpindex x nvm words = ',&
1268            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm 
1269       STOP 'sechiba_init'
1270    END IF
1271
1272    ALLOCATE (rveget(kjpindex,nvm),stat=ier)
1273    IF (ier.NE.0) THEN
1274       WRITE (numout,*) ' error in rveget allocation. We stop. We need kjpindex x nvm words = ',&
1275            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1276       STOP 'sechiba_init'
1277    END IF
1278
1279    !
1280    !  1.5  Get the indexing table for the vegetation fields. In SECHIBA we work on reduced grids but to store in the
1281    !          full 3D filed vegetation variable we need another index table : indexveg, indexsoil, indexnobio and
1282    !          indexgrnd
1283    !
1284    DO ji = 1, kjpindex
1285       !
1286       DO jv = 1, nvm
1287          indexveg((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1288       ENDDO
1289       !     
1290       DO jv = 1, nstm
1291          indexsoil((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1292       ENDDO
1293       !     
1294       DO jv = 1, nnobio
1295          indexnobio((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1296       ENDDO
1297       !
1298       DO jv = 1, ngrnd
1299          indexgrnd((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1300       ENDDO
1301       !
1302       DO jv = 1, nslm
1303          indexlayer((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1304       ENDDO
1305       !
1306       DO jv = 1, 2
1307          indexalb((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1308       ENDDO
1309       !
1310    ENDDO
1311
1312    !
1313    ! 2. restart value
1314    !
1315    ! open restart input file if needed
1316    ! and read data from restart input file
1317    !
1318
1319    IF (ldrestart_read) THEN
1320
1321       IF (long_print) WRITE (numout,*) ' we have to read a restart file for SECHIBA variables'
1322       !
1323       !   Read the default value that will be put into variable which are not in the restart file
1324       !
1325       CALL ioget_expval(val_exp)
1326
1327    ENDIF
1328
1329    !
1330    river_routing = control_in%river_routing
1331    hydrol_cwrr = control_in%hydrol_cwrr
1332    !
1333    ! 4. run control: store flags in a common variable
1334    !
1335
1336    control%ok_co2 = control_in%ok_co2
1337    control%ok_sechiba = control_in%ok_sechiba
1338    control%ok_stomate = control_in%ok_stomate
1339    control%ok_dgvm = control_in%ok_dgvm
1340    control%ok_pheno = control_in%ok_pheno
1341    control%stomate_watchout = control_in%stomate_watchout
1342
1343    IF (long_print) WRITE (numout,*) ' sechiba_init done '
1344
1345  END SUBROUTINE sechiba_init
1346  !
1347  !------------------------------------------------------------------
1348  !
1349  SUBROUTINE sechiba_clear (forcing_name,cforcing_name)
1350
1351    CHARACTER(LEN=100), INTENT(in)           :: forcing_name
1352    CHARACTER(LEN=100), INTENT(in)           :: cforcing_name
1353
1354    !
1355    ! initialisation
1356    !
1357    l_first_sechiba=.TRUE.
1358
1359    ! 1. Deallocate all dynamic variables
1360
1361    IF ( ALLOCATED (indexveg)) DEALLOCATE (indexveg)
1362    IF ( ALLOCATED (indexsoil)) DEALLOCATE (indexsoil)
1363    IF ( ALLOCATED (indexnobio)) DEALLOCATE (indexnobio)
1364    IF ( ALLOCATED (indexgrnd)) DEALLOCATE (indexgrnd)
1365    IF ( ALLOCATED (indexlayer)) DEALLOCATE (indexlayer)
1366    IF ( ALLOCATED (indexalb)) DEALLOCATE (indexalb)
1367    IF ( ALLOCATED (snow)) DEALLOCATE (snow)
1368    IF ( ALLOCATED (snow_age)) DEALLOCATE (snow_age)
1369    IF ( ALLOCATED (drysoil_frac)) DEALLOCATE (drysoil_frac)
1370    IF ( ALLOCATED (rsol)) DEALLOCATE (rsol)
1371    IF ( ALLOCATED (evap_bare_lim)) DEALLOCATE (evap_bare_lim)
1372    IF ( ALLOCATED (evapot)) DEALLOCATE (evapot)
1373    IF ( ALLOCATED (evapot_corr)) DEALLOCATE (evapot_corr)
1374    IF ( ALLOCATED (humrel)) DEALLOCATE (humrel)
1375    IF ( ALLOCATED (vegstress)) DEALLOCATE (vegstress)
1376    IF ( ALLOCATED (soiltype)) DEALLOCATE (soiltype)
1377    IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1)
1378    IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4)
1379    IF ( ALLOCATED (soilcap)) DEALLOCATE (soilcap)
1380    IF ( ALLOCATED (soilflx)) DEALLOCATE (soilflx)
1381    IF ( ALLOCATED (temp_sol)) DEALLOCATE (temp_sol)
1382    IF ( ALLOCATED (qsurf)) DEALLOCATE (qsurf)
1383    IF ( ALLOCATED (qsintveg)) DEALLOCATE (qsintveg)
1384    IF ( ALLOCATED (vbeta2))  DEALLOCATE (vbeta2)
1385    IF ( ALLOCATED (vbeta3)) DEALLOCATE (vbeta3)
1386    IF ( ALLOCATED (vbetaco2)) DEALLOCATE (vbetaco2)
1387    IF ( ALLOCATED (cimean)) DEALLOCATE (cimean)
1388    IF ( ALLOCATED (gpp)) DEALLOCATE (gpp)
1389    IF ( ALLOCATED (veget)) DEALLOCATE (veget)
1390    IF ( ALLOCATED (veget_max)) DEALLOCATE (veget_max)
1391    IF ( ALLOCATED (lai)) DEALLOCATE (lai)
1392    IF ( ALLOCATED (height)) DEALLOCATE (height)
1393    IF ( ALLOCATED (roughheight)) DEALLOCATE (roughheight)
1394    IF ( ALLOCATED (frac_nobio)) DEALLOCATE (frac_nobio)
1395    IF ( ALLOCATED (snow_nobio)) DEALLOCATE (snow_nobio)
1396    IF ( ALLOCATED (snow_nobio_age)) DEALLOCATE (snow_nobio_age)
1397    IF ( ALLOCATED (assim_param)) DEALLOCATE (assim_param)
1398    IF ( ALLOCATED (vevapsno)) DEALLOCATE (vevapsno)
1399    IF ( ALLOCATED (vevapnu)) DEALLOCATE (vevapnu)
1400    IF ( ALLOCATED (t2mdiag)) DEALLOCATE (t2mdiag)
1401    IF ( ALLOCATED (totfrac_nobio)) DEALLOCATE (totfrac_nobio)
1402    IF ( ALLOCATED (runoff)) DEALLOCATE (runoff)
1403    IF ( ALLOCATED (drainage)) DEALLOCATE (drainage)
1404    IF ( ALLOCATED (returnflow)) DEALLOCATE (returnflow)
1405    IF ( ALLOCATED (irrigation)) DEALLOCATE (irrigation)
1406    IF ( ALLOCATED (tot_melt)) DEALLOCATE (tot_melt)
1407    IF ( ALLOCATED (valpha)) DEALLOCATE (valpha)
1408    IF ( ALLOCATED (vbeta)) DEALLOCATE (vbeta)
1409    IF ( ALLOCATED (fusion)) DEALLOCATE (fusion)
1410    IF ( ALLOCATED (rau)) DEALLOCATE (rau)
1411    IF ( ALLOCATED (deadleaf_cover)) DEALLOCATE (deadleaf_cover)
1412    IF ( ALLOCATED (stempdiag)) DEALLOCATE (stempdiag)
1413    IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
1414    IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
1415    IF ( ALLOCATED (litterhumdiag)) DEALLOCATE (litterhumdiag)
1416    IF ( ALLOCATED (vevapwet)) DEALLOCATE (vevapwet)
1417    IF ( ALLOCATED (transpir)) DEALLOCATE (transpir)
1418    IF ( ALLOCATED (qsintmax)) DEALLOCATE (qsintmax)
1419    IF ( ALLOCATED (rveget)) DEALLOCATE (rveget)
1420    ! 2. clear all modules
1421    CALL slowproc_clear 
1422    CALL diffuco_clear 
1423    CALL enerbil_clear 
1424    IF ( hydrol_cwrr ) THEN
1425       CALL hydrol_clear 
1426    ELSE
1427       CALL hydrolc_clear 
1428    ENDIF
1429    CALL condveg_clear 
1430    CALL thermosoil_clear
1431    CALL routing_clear
1432    !3. give name to next block
1433    stomate_forcing_name=forcing_name
1434    stomate_Cforcing_name=Cforcing_name
1435
1436  END SUBROUTINE sechiba_clear
1437
1438  !! SECHIBA's variables initialisation
1439  !! called every time step
1440  !!
1441  SUBROUTINE sechiba_var_init (kjpindex, rau, pb, temp_air) 
1442
1443    ! interface description
1444    ! input scalar
1445    INTEGER(i_std), INTENT (in)                    :: kjpindex           !! Domain dimension
1446    ! input fields
1447    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: pb                 !! Lowest level pressure
1448    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: temp_air           !! Air temperature
1449    ! output fields
1450    REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rau                !! Density
1451
1452    ! local declaration
1453    INTEGER(i_std)                                :: ji
1454
1455    !
1456    ! initialisation
1457    !
1458
1459    !
1460    ! 1. calcul of rau: air density
1461    !
1462
1463    DO ji = 1,kjpindex
1464       rau(ji) = pa_par_hpa * pb(ji) / (cte_molr*temp_air(ji))
1465    END DO
1466
1467    IF (long_print) WRITE (numout,*) ' sechiba_var_init done '
1468
1469  END SUBROUTINE sechiba_var_init
1470
1471  !!
1472  !! Swap new fields to previous fields
1473  !!
1474  SUBROUTINE sechiba_end (kjpindex, dtradia, temp_sol, temp_sol_new)
1475
1476    ! interface description
1477    ! input scalar
1478    INTEGER(i_std), INTENT (in)                       :: kjpindex           !! Domain dimension
1479    REAL(r_std),INTENT (in)                           :: dtradia            !! Time step in seconds
1480    ! input fields
1481    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: temp_sol_new       !! New soil temperature
1482    ! output fields
1483    REAL(r_std),DIMENSION (kjpindex), INTENT (out)    :: temp_sol           !! Soil temperature
1484
1485    !
1486    ! swap
1487    !
1488    temp_sol(:) = temp_sol_new(:)
1489
1490    IF (long_print) WRITE (numout,*) ' sechiba_end done '
1491
1492  END SUBROUTINE sechiba_end
1493
1494END MODULE sechiba
Note: See TracBrowser for help on using the repository browser.