source: branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/sechiba.f90 @ 143

Last change on this file since 143 was 142, checked in by didier.solyga, 13 years ago

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