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

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

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