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