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

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

Externalized version merged with the trunk
File size: 69.1 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: 45 $, $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $
7!!
8!< $HeadURL: http://forge.ipsl.jussieu.fr/orchidee/svn/trunk/ORCHIDEE/src_sechiba/sechiba.f90 $
9!< $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $
10!< $Author: mmaipsl $
11!< $Revision: 45 $
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, &
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
252    REAL(r_std), ALLOCATABLE, DIMENSION (:)                  :: runoff1,drainage1, soilcap1,soilflx1
253    REAL(r_std), ALLOCATABLE, DIMENSION (:,:)                :: shumdiag1
254
255    REAL(r_std), DIMENSION(kjpindex)                   :: histvar          !! computations for history files
256
257
258    ! 15/03/2011 DS Add for externalisation
259    REAL(r_std), DIMENSION(kjpindex) :: sum_treefrac, sum_grassfrac, sum_cropfrac
260    INTEGER(i_std) :: jv
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)/one_day
639       CALL histwrite(hist_id, 'evspsblveg', kjit, histvar, kjpindex, index)
640
641       histvar(:)=(vevapnu(:)+vevapsno(:))/one_day
642       CALL histwrite(hist_id, 'evspsblsoi', kjit, histvar, kjpindex, index)
643
644       histvar(:)=SUM(transpir(:,:),dim=2)/one_day
645       CALL histwrite(hist_id, 'tran', kjit, histvar, kjpindex, index)
646
647!$$ 25/10/10 Modif DS & NViovy
648!!$
649       histvar(:)= sum_treefrac(:)*100*contfrac(:)
650       CALL histwrite(hist_id, 'treeFrac', kjit, histvar, kjpindex, index) 
651
652       histvar(:)= sum_grassfrac(:)*100*contfrac(:)
653       CALL histwrite(hist_id, 'grassFrac', kjit, histvar, kjpindex, index) 
654
655       histvar(:)= sum_cropfrac(:)*100*contfrac(:)
656       CALL histwrite(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)
657
658       histvar(:)=veget_max(:,1)*100*contfrac(:)
659       CALL histwrite(hist_id, 'baresoilFrac', kjit, histvar, kjpindex, index)
660
661       histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)*100*contfrac(:)
662       CALL histwrite(hist_id, 'residualFrac', kjit, histvar, kjpindex, index)
663
664    ELSE
665       CALL histwrite(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
666       CALL histwrite(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
667       CALL histwrite(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
668       CALL histwrite(hist_id, 'Qf', kjit, fusion, kjpindex, index)
669       CALL histwrite(hist_id, 'SWE', kjit, snow, kjpindex, index)
670       CALL histwrite(hist_id, 'ESoil', kjit, vevapnu, kjpindex, index)
671       CALL histwrite(hist_id, 'TVeg', kjit, transpir, kjpindex*nvm, indexveg)
672       CALL histwrite(hist_id, 'ECanop', kjit, vevapwet, kjpindex*nvm, indexveg)
673       CALL histwrite(hist_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
674       CALL histwrite(hist_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
675       IF ( control%ok_co2 ) THEN
676          CALL histwrite(hist_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
677       ENDIF
678       IF ( control%ok_stomate ) THEN
679             CALL histwrite(hist_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)   
680       ENDIF
681    ENDIF
682    !
683    IF ( hist2_id > 0 ) THEN
684       IF ( .NOT. almaoutput ) THEN
685          CALL histwrite(hist2_id, 'tsol_rad', kjit, tsol_rad, kjpindex, index)
686          CALL histwrite(hist2_id, 'qsurf', kjit, qsurf, kjpindex, index)
687          CALL histwrite(hist2_id, 'albedo', kjit, albedo, kjpindex*2, indexalb)
688          CALL histwrite(hist2_id, 'emis', kjit, emis, kjpindex, index)
689          !
690          CALL histwrite(hist2_id, 'beta', kjit, vbeta, kjpindex, index)
691          CALL histwrite(hist2_id, 'z0', kjit, z0, kjpindex, index)
692          CALL histwrite(hist2_id, 'roughheight', kjit, roughheight, kjpindex, index)
693          CALL histwrite(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
694          CALL histwrite(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
695          CALL histwrite(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
696          CALL histwrite(hist2_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
697          CALL histwrite(hist2_id, 'subli', kjit, vevapsno, kjpindex, index)
698          CALL histwrite(hist2_id, 'vevapnu', kjit, vevapnu, kjpindex, index)
699          CALL histwrite(hist2_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
700          CALL histwrite(hist2_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
701          CALL histwrite(hist2_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
702          CALL histwrite(hist2_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
703          CALL histwrite(hist2_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
704          CALL histwrite(hist2_id, 'vbeta4', kjit, vbeta4, kjpindex, index)   
705          CALL histwrite(hist2_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
706          CALL histwrite(hist2_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
707          CALL histwrite(hist2_id, 'rsol', kjit, rsol, kjpindex, index)
708          CALL histwrite(hist2_id, 'snow', kjit, snow, kjpindex, index)
709          CALL histwrite(hist2_id, 'snowage', kjit, snow_age, kjpindex, index)
710          CALL histwrite(hist2_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
711          CALL histwrite(hist2_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
712          CALL histwrite(hist2_id, 'soiltype',  kjit, soiltype, kjpindex*nstm, indexsoil)
713          IF ( control%ok_co2 ) THEN
714             CALL histwrite(hist2_id, 'vbetaco2', kjit, vbetaco2, kjpindex*nvm, indexveg)   
715             CALL histwrite(hist2_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
716             CALL histwrite(hist2_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)   
717          ENDIF
718          IF ( control%ok_stomate ) THEN
719             CALL histwrite(hist2_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)   
720          ENDIF
721       ELSE
722          CALL histwrite(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
723          CALL histwrite(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
724          CALL histwrite(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
725          CALL histwrite(hist2_id, 'Qf', kjit, fusion, kjpindex, index)
726          CALL histwrite(hist2_id, 'SWE', kjit, snow, kjpindex, index)
727          CALL histwrite(hist2_id, 'ESoil', kjit, vevapnu, kjpindex, index)
728          CALL histwrite(hist2_id, 'TVeg', kjit, transpir, kjpindex*nvm, indexveg)
729          CALL histwrite(hist2_id, 'ECanop', kjit, vevapwet, kjpindex*nvm, indexveg)
730          CALL histwrite(hist2_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
731          CALL histwrite(hist2_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
732          IF ( control%ok_co2 ) THEN
733             CALL histwrite(hist2_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
734          ENDIF
735          IF ( control%ok_stomate ) THEN
736             CALL histwrite(hist2_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)   
737          ENDIF
738       ENDIF
739    ENDIF
740    !
741    ! prepares restart file for the next simulation from SECHIBA and from other modules
742    !
743    IF (ldrestart_write) THEN
744
745       IF (long_print) WRITE (numout,*) ' we have to write a restart file '
746
747       !
748       ! call diffuco_main to write restart files
749       !
750
751       CALL diffuco_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, u, v, &
752! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
753!           & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, pb ,  &
754            & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, q2m, t2m, pb ,  &
755            & rsol, evap_bare_lim, evapot, snow, frac_nobio, snow_nobio, totfrac_nobio, &
756            & swnet, swdown, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
757            & vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, rveget, cimean, rest_id, hist_id, hist2_id)
758
759       !
760       ! call energy bilan to write restart files
761       !
762       CALL enerbil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, &
763            & index, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef,&
764            & qair, peqAcoef, peqBcoef, pb, rau, vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, &
765            & cimean, ccanopy, emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
766            & vevapp, transpir, gpp, vevapnu, vevapwet, vevapsno, t2mdiag, temp_sol, tsol_rad, &
767            & temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id) 
768
769       !
770       ! call hydrologie to write restart files
771       !
772       IF ( .NOT. hydrol_cwrr ) THEN
773          CALL hydrolc_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, &
774               & temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
775               & qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
776               & tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, &
777               & humrel, vegstress, rsol, drysoil_frac, evapot, evapot_corr, shumdiag, litterhumdiag, soilcap, &
778               & rest_id, hist_id, hist2_id)
779          evap_bare_lim(:) = -un
780       ELSE
781          CALL hydrol_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, indexsoil, indexlayer, &
782               & temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
783               & qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
784               & tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
785               & vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, &
786               & shumdiag, litterhumdiag, soilcap, soiltype, rest_id, hist_id, hist2_id)
787          rsol(:) = -un
788       ENDIF
789       !
790       ! call condveg to write restart files
791       !
792       CALL condveg_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, &
793            & lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
794            & zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
795            & drysoil_frac, height, deadleaf_cover, emis, albedo, z0, roughheight, soiltype, rest_id, hist_id, hist2_id)
796
797       !
798       ! call Soil Thermodynamic to write restart files
799       !
800       CALL thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexgrnd, &
801            & temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id) 
802
803       !
804       ! If we chose to route the water then we call the module. Else variables
805       ! are set to zero.
806       !
807       !
808       IF ( river_routing .AND. nbp_glo .GT. 1) THEN
809          CALL routing_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, &
810               & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, runoff, &
811               & drainage, evapot_corr, precip_rain, humrel, &
812               & stempdiag, returnflow, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
813       ELSE
814          riverflow(:) = zero
815          coastalflow(:) = zero
816          returnflow(:) = zero
817          irrigation(:) = zero
818       ENDIF
819
820       !
821       ! call slowproc_main to write restart files
822       !
823
824       CALL slowproc_main (kjit, kjpij, kjpindex, dtradia, date0, &
825            ldrestart_read, ldrestart_write, control%ok_co2, control%ok_stomate, &
826            index, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
827            t2mdiag, t2mdiag, temp_sol, stempdiag, &
828            vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
829            deadleaf_cover, &
830            assim_param, &
831            lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
832            rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
833            co2_flux)
834
835
836       var_name= 'shumdiag' 
837       CALL restput_p(rest_id, var_name, nbp_glo,   nbdl, 1, kjit,  shumdiag, 'scatter',  nbp_glo, index_g)
838
839       var_name= 'runoff' 
840       CALL restput_p(rest_id, var_name, nbp_glo,   1, 1, kjit,  runoff, 'scatter',  nbp_glo, index_g)
841       var_name= 'drainage' 
842       CALL restput_p(rest_id, var_name, nbp_glo,   1, 1, kjit,  drainage, 'scatter',  nbp_glo, index_g)
843    END IF
844
845    IF (long_print) WRITE (numout,*) ' sechiba_main done '
846
847  END SUBROUTINE sechiba_main
848
849  !! Initialisation for SECHIBA processes
850  !! - does dynamic allocation for local arrays
851  !! - reads _restart_ file or set initial values to a raisonable value
852  !! - reads initial map
853  !!
854  SUBROUTINE sechiba_init (kjit, ldrestart_read, kjpij, kjpindex, index, rest_id, control_in, lalo)
855
856    ! interface description
857    ! input scalar
858    INTEGER(i_std), INTENT (in)                         :: kjit               !! Time step number
859    LOGICAL,INTENT (in)                                :: ldrestart_read     !! Logical for restart file to read
860    INTEGER(i_std), INTENT (in)                         :: kjpij              !! Size of full domaine
861    INTEGER(i_std), INTENT (in)                         :: kjpindex           !! Domain size
862    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: index              !! Indeces of the points on the map
863    INTEGER(i_std), INTENT (in)                         :: rest_id            !! _Restart_ file identifier
864    TYPE(control_type), INTENT(in)                     :: control_in         !! Flags that (de)activate parts of the model
865    ! input fields
866    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)     :: lalo               !! Geographical coordinates
867    ! output scalar
868    ! output fields
869
870    ! local declaration
871    INTEGER(i_std)                                :: ier, ji, jv
872    !
873    ! initialisation
874    !
875    IF (l_first_sechiba) THEN
876       l_first_sechiba=.FALSE.
877    ELSE
878       WRITE (numout,*) ' l_first_sechiba false . we stop '
879       STOP 'sechiba_init'
880    ENDIF
881
882    ! 1. make dynamic allocation with good dimension
883
884    ! 1.0 The 3D vegetation indexation table
885
886    ALLOCATE (indexveg(kjpindex*nvm),stat=ier)
887    IF (ier.NE.0) THEN
888       WRITE (numout,*) ' error in indexveg allocation. We stop. We need kjpindex words = ',kjpindex*nvm
889       STOP 'sechiba_init'
890    END IF
891
892    ALLOCATE (indexsoil(kjpindex*nstm),stat=ier)
893    IF (ier.NE.0) THEN
894       WRITE (numout,*) ' error in indexsoil allocation. We stop. We need kjpindex words = ',kjpindex*nstm
895       STOP 'sechiba_init'
896    END IF
897
898    ALLOCATE (indexnobio(kjpindex*nnobio),stat=ier)
899    IF (ier.NE.0) THEN
900       WRITE (numout,*) ' error in indexnobio allocation. We stop. We need kjpindex words = ',kjpindex*nnobio
901       STOP 'sechiba_init'
902    END IF
903
904    ALLOCATE (indexgrnd(kjpindex*ngrnd),stat=ier)
905    IF (ier.NE.0) THEN
906       WRITE (numout,*) ' error in indexgrnd allocation. We stop. We need kjpindex words = ',kjpindex*ngrnd
907       STOP 'sechiba_init'
908    END IF
909
910    ALLOCATE (indexlayer(kjpindex*nslm),stat=ier)
911    IF (ier.NE.0) THEN
912       WRITE (numout,*) ' error in indexlayer allocation. We stop. We need kjpindex words = ',kjpindex*nslm
913       STOP 'sechiba_init'
914    END IF
915
916    ALLOCATE (indexalb(kjpindex*2),stat=ier)
917    IF (ier.NE.0) THEN
918       WRITE (numout,*) ' error in indexalb allocation. We stop. We need kjpindex words = ',kjpindex*2
919       STOP 'sechiba_init'
920    END IF
921
922    ! 1.1 one dimension array allocation with restartable value
923
924    IF (long_print) WRITE (numout,*) 'Allocation of 1D variables. We need for each kjpindex words = ',kjpindex
925    ALLOCATE (snow(kjpindex),stat=ier)
926    IF (ier.NE.0) THEN
927       WRITE (numout,*) ' error in snow allocation. We stop. We need kjpindex words = ',kjpindex
928       STOP 'sechiba_init'
929    END IF
930    snow(:) = undef_sechiba
931
932    ALLOCATE (snow_age(kjpindex),stat=ier)
933    IF (ier.NE.0) THEN
934       WRITE (numout,*) ' error in snow_age allocation. We stop. We need kjpindex words = ',kjpindex
935       STOP 'sechiba_init'
936    END IF
937    snow_age(:) = undef_sechiba
938
939    ALLOCATE (drysoil_frac(kjpindex),stat=ier)
940    IF (ier.NE.0) THEN
941       WRITE (numout,*) ' error in drysoil_frac allocation. We stop. We need kjpindex words = ',kjpindex
942       STOP 'sechiba_init'
943    END IF
944    drysoil_frac(:) = zero
945
946    ALLOCATE (rsol(kjpindex),stat=ier)
947    IF (ier.NE.0) THEN
948       WRITE (numout,*) ' error in rsol allocation. We stop. We need kjpindex words = ',kjpindex
949       STOP 'sechiba_init'
950    END IF
951
952    ALLOCATE (evap_bare_lim(kjpindex),stat=ier)
953    IF (ier.NE.0) THEN
954       WRITE (numout,*) ' error in evap_bare_lim allocation. We stop. We need kjpindex words = ',kjpindex
955       STOP 'sechiba_init'
956    END IF
957
958    ALLOCATE (evapot(kjpindex),stat=ier)
959    IF (ier.NE.0) THEN
960       WRITE (numout,*) ' error in evapot allocation. We stop. We need kjpindex words = ',kjpindex
961       STOP 'sechiba_init'
962    END IF
963    evapot(:) = undef_sechiba
964
965    ALLOCATE (evapot_corr(kjpindex),stat=ier)
966    IF (ier.NE.0) THEN
967       WRITE (numout,*) ' error in evapot_corr allocation. We stop. We need kjpindex words = ',kjpindex
968       STOP 'sechiba_init'
969    END IF
970
971    ALLOCATE (humrel(kjpindex,nvm),stat=ier)
972    IF (ier.NE.0) THEN
973       WRITE (numout,*) ' error in humrel allocation. We stop. we need kjpindex words = ',kjpindex
974       STOP 'sechiba_init'
975    END IF
976    humrel(:,:) = undef_sechiba
977
978    ALLOCATE (vegstress(kjpindex,nvm),stat=ier)
979    IF (ier.NE.0) THEN
980       WRITE (numout,*) ' error in vegstress allocation. We stop. we need kjpindex words = ',kjpindex
981       STOP 'sechiba_init'
982    END IF
983    vegstress(:,:) = undef_sechiba
984
985    ALLOCATE (soiltype(kjpindex,nstm),stat=ier)
986    IF (ier.NE.0) THEN
987       WRITE (numout,*) ' error in soiltype allocation. We stop. we need kjpindex words = ',kjpindex
988       STOP 'sechiba_init'
989    END IF
990    soiltype(:,:)=undef_sechiba
991
992    ALLOCATE (vbeta1(kjpindex),stat=ier)
993    IF (ier.NE.0) THEN
994       WRITE (numout,*) ' error in vbeta1 allocation. We stop. We need kjpindex words = ',kjpindex
995       STOP 'sechiba_init'
996    END IF
997
998    ALLOCATE (vbeta4(kjpindex),stat=ier)
999    IF (ier.NE.0) THEN
1000       WRITE (numout,*) ' error in vbeta4 allocation. We stop. We need kjpindex words = ',kjpindex
1001       STOP 'sechiba_init'
1002    END IF
1003
1004    ALLOCATE (soilcap(kjpindex),stat=ier)
1005    IF (ier.NE.0) THEN
1006       WRITE (numout,*) ' error in soilcap allocation. We stop. We need kjpindex words = ',kjpindex
1007       STOP 'sechiba_init'
1008    END IF
1009
1010    ALLOCATE (soilflx(kjpindex),stat=ier)
1011    IF (ier.NE.0) THEN
1012       WRITE (numout,*) ' error in soilflx allocation. We stop. We need kjpindex words = ',kjpindex
1013       STOP 'sechiba_init'
1014    END IF
1015
1016    ALLOCATE (temp_sol(kjpindex),stat=ier)
1017    IF (ier.NE.0) THEN
1018       WRITE (numout,*) ' error in temp_sol allocation. We stop. We need kjpindex words = ',kjpindex
1019       STOP 'sechiba_init'
1020    END IF
1021    temp_sol(:) = undef_sechiba
1022
1023    ALLOCATE (qsurf(kjpindex),stat=ier)
1024    IF (ier.NE.0) THEN
1025       WRITE (numout,*) ' error in qsurf allocation. We stop. We need kjpindex words = ',kjpindex
1026       STOP 'sechiba_init'
1027    END IF
1028    qsurf(:) = undef_sechiba
1029
1030    ! 1.2 two dimensions array allocation with restartable value
1031
1032    ALLOCATE (qsintveg(kjpindex,nvm),stat=ier)
1033    IF (ier.NE.0) THEN
1034       WRITE (numout,*) ' error in qsintveg allocation. We stop. We need kjpindex x nvm words = ',&
1035            & kjpindex,' x ' ,nvm,' = ',kjpindex*nvm 
1036       STOP 'sechiba_init'
1037    END IF
1038    qsintveg(:,:) = undef_sechiba
1039
1040    ALLOCATE (vbeta2(kjpindex,nvm),stat=ier)
1041    IF (ier.NE.0) THEN
1042       WRITE (numout,*) ' error in vbeta2 allocation. We stop. We need kjpindex x nvm words = ',&
1043            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1044       STOP 'sechiba_init'
1045    END IF
1046
1047    ALLOCATE (vbeta3(kjpindex,nvm),stat=ier)
1048    IF (ier.NE.0) THEN
1049       WRITE (numout,*) ' error in vbeta3 allocation. We stop.We need kjpindex x nvm words = ',&
1050            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1051       STOP 'sechiba_init'
1052    END IF
1053
1054    ALLOCATE (vbetaco2(kjpindex,nvm),stat=ier)
1055    IF (ier.NE.0) THEN
1056       WRITE (numout,*) ' error in vbetaco2 allocation. We stop.We need kjpindex x nvm words = ',&
1057            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1058       STOP 'sechiba_init'
1059    END IF
1060
1061    ALLOCATE (cimean(kjpindex,nvm),stat=ier)
1062    IF (ier.NE.0) THEN
1063       WRITE (numout,*) ' error in cimean allocation. We stop.We need kjpindex x nvm words = ',&
1064            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1065       STOP 'sechiba_init'
1066    END IF
1067
1068    ALLOCATE (gpp(kjpindex,nvm),stat=ier)
1069    IF (ier.NE.0) THEN
1070       WRITE (numout,*) ' error in gpp allocation. We stop.We need kjpindex x nvm words = ',&
1071            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1072       STOP 'sechiba_init'
1073    END IF
1074    gpp(:,:) = undef_sechiba
1075
1076    ALLOCATE (veget(kjpindex,nvm),stat=ier)
1077    IF (ier.NE.0) THEN
1078       WRITE (numout,*) ' error in veget allocation. We stop. We need kjpindex x nvm words = ',&
1079            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1080       STOP 'sechiba_init'
1081    END IF
1082    veget(:,:)=undef_sechiba
1083
1084    ALLOCATE (veget_max(kjpindex,nvm),stat=ier)
1085    IF (ier.NE.0) THEN
1086       WRITE (numout,*) ' error in veget_max allocation. We stop. We need kjpindex x nvm words = ',&
1087            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1088       STOP 'sechiba_init'
1089    END IF
1090    veget_max(:,:)=undef_sechiba
1091
1092    ALLOCATE (lai(kjpindex,nvm),stat=ier)
1093    IF (ier.NE.0) THEN
1094       WRITE (numout,*) ' error in lai allocation. We stop. We need kjpindex x nvm words = ',&
1095            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1096       STOP 'sechiba_init'
1097    END IF
1098    lai(:,:)=undef_sechiba
1099
1100    ALLOCATE (height(kjpindex,nvm),stat=ier)
1101    IF (ier.NE.0) THEN
1102       WRITE (numout,*) ' error in height allocation. We stop. We need kjpindex x nvm words = ',&
1103            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1104       STOP 'sechiba_init'
1105    END IF
1106    height(:,:)=undef_sechiba
1107
1108    ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
1109    IF (ier.NE.0) THEN
1110       WRITE (numout,*) ' error in frac_nobio allocation. We stop. We need kjpindex words = ',kjpindex
1111       STOP 'sechiba_init'
1112    END IF
1113    frac_nobio(:,:) = undef_sechiba
1114
1115    ALLOCATE (albedo(kjpindex,2),stat=ier)
1116    IF (ier.NE.0) THEN
1117       WRITE (numout,*) ' error in albedo allocation. We stop. We need kjpindex words = ',kjpindex
1118       STOP 'sechiba_init'
1119    END IF
1120
1121    ALLOCATE (snow_nobio(kjpindex,nnobio),stat=ier)
1122    IF (ier.NE.0) THEN
1123       WRITE (numout,*) ' error in snow_nobio allocation. We stop. We need kjpindex words = ',kjpindex
1124       STOP 'sechiba_init'
1125    END IF
1126    snow_nobio(:,:) = undef_sechiba
1127
1128    ALLOCATE (snow_nobio_age(kjpindex,nnobio),stat=ier)
1129    IF (ier.NE.0) THEN
1130       WRITE (numout,*) ' error in snow_nobio_age allocation. We stop. We need kjpindex words = ',kjpindex
1131       STOP 'sechiba_init'
1132    END IF
1133    snow_nobio_age(:,:) = undef_sechiba
1134
1135    ALLOCATE (assim_param(kjpindex,nvm,npco2),stat=ier)
1136    IF (ier.NE.0) THEN
1137       WRITE (numout,*) ' error in assim_param allocation. We stop. We need kjpindex x nvm x npco2 words = ',&
1138            & kjpindex,' x ' ,nvm,' x ',npco2, ' = ',kjpindex*nvm*npco2
1139       STOP 'sechiba_init'
1140    END IF
1141
1142    ! 1.3 one dimension array allocation
1143
1144    ALLOCATE (vevapsno(kjpindex),stat=ier)
1145    IF (ier.NE.0) THEN
1146       WRITE (numout,*) ' error in vevapsno allocation. We stop. We need kjpindex words = ',kjpindex
1147       STOP 'sechiba_init'
1148    END IF
1149
1150    ALLOCATE (vevapnu(kjpindex),stat=ier)
1151    IF (ier.NE.0) THEN
1152       WRITE (numout,*) ' error in vevapnu allocation. We stop. We need kjpindex words = ',kjpindex
1153       STOP 'sechiba_init'
1154    END IF
1155
1156    ALLOCATE (t2mdiag(kjpindex),stat=ier)
1157    IF (ier.NE.0) THEN
1158       WRITE (numout,*) ' error in t2mdiag allocation. We stop. We need kjpindex words = ',kjpindex
1159       STOP 'sechiba_init'
1160    END IF
1161
1162    ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
1163    IF (ier.NE.0) THEN
1164       WRITE (numout,*) ' error in totfrac_nobio allocation. We stop. We need kjpindex words = ',kjpindex
1165       STOP 'sechiba_init'
1166    END IF
1167
1168    ALLOCATE (runoff(kjpindex),stat=ier)
1169    IF (ier.NE.0) THEN
1170       WRITE (numout,*) ' error in runoff allocation. We stop. We need kjpindex words = ',kjpindex
1171       STOP 'sechiba_init'
1172    END IF
1173
1174    ALLOCATE (drainage(kjpindex),stat=ier)
1175    IF (ier.NE.0) THEN
1176       WRITE (numout,*) ' error in drainage allocation. We stop. We need kjpindex words = ',kjpindex
1177       STOP 'sechiba_init'
1178    END IF
1179
1180    ALLOCATE (returnflow(kjpindex),stat=ier)
1181    IF (ier.NE.0) THEN
1182       WRITE (numout,*) ' error in returnflow allocation. We stop. We need kjpindex words = ',kjpindex
1183       STOP 'sechiba_init'
1184    END IF
1185    returnflow(:) = zero
1186
1187    ALLOCATE (irrigation(kjpindex),stat=ier)
1188    IF (ier.NE.0) THEN
1189       WRITE (numout,*) ' error in irrigation allocation. We stop. We need kjpindex words = ',kjpindex
1190       STOP 'sechiba_init'
1191    END IF
1192    irrigation(:) = zero
1193
1194    ALLOCATE (z0(kjpindex),stat=ier)
1195    IF (ier.NE.0) THEN
1196       WRITE (numout,*) ' error in z0 allocation. We stop. We need kjpindex words = ',kjpindex
1197       STOP 'sechiba_init'
1198    END IF
1199
1200    ALLOCATE (roughheight(kjpindex),stat=ier)
1201    IF (ier.NE.0) THEN
1202       WRITE (numout,*) ' error in roughheight allocation. We stop. We need kjpindex words = ',kjpindex
1203       STOP 'sechiba_init'
1204    END IF
1205
1206    ALLOCATE (emis(kjpindex),stat=ier)
1207    IF (ier.NE.0) THEN
1208       WRITE (numout,*) ' error in emis allocation. We stop. We need kjpindex words = ',kjpindex
1209       STOP 'sechiba_init'
1210    END IF
1211
1212    ALLOCATE (tot_melt(kjpindex),stat=ier)
1213    IF (ier.NE.0) THEN
1214       WRITE (numout,*) ' error in tot_melt allocation. We stop. We need kjpindex words = ',kjpindex
1215       STOP 'sechiba_init'
1216    END IF
1217
1218    ALLOCATE (valpha(kjpindex),stat=ier)
1219    IF (ier.NE.0) THEN
1220       WRITE (numout,*) ' error in valpha allocation. We stop. We need kjpindex words = ',kjpindex
1221       STOP 'sechiba_init'
1222    END IF
1223
1224    ALLOCATE (vbeta(kjpindex),stat=ier)
1225    IF (ier.NE.0) THEN
1226       WRITE (numout,*) ' error in vbeta allocation. We stop. We need kjpindex words = ',kjpindex
1227       STOP 'sechiba_init'
1228    END IF
1229
1230    ALLOCATE (fusion(kjpindex),stat=ier)
1231    IF (ier.NE.0) THEN
1232       WRITE (numout,*) ' error in fusion allocation. We stop. We need kjpindex words = ',kjpindex
1233       STOP 'sechiba_init'
1234    END IF
1235
1236    ALLOCATE (rau(kjpindex),stat=ier)
1237    IF (ier.NE.0) THEN
1238       WRITE (numout,*) ' error in rau allocation. We stop. We need kjpindex words = ',kjpindex
1239       STOP 'sechiba_init'
1240    END IF
1241
1242    ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
1243    IF (ier.NE.0) THEN
1244       WRITE (numout,*) ' error in deadleaf_cover allocation. We stop. We need kjpindex words = ',kjpindex
1245       STOP 'sechiba_init'
1246    END IF
1247
1248    ALLOCATE (stempdiag(kjpindex, nbdl),stat=ier)
1249    IF (ier.NE.0) THEN
1250       WRITE (numout,*) ' error in stempdiag allocation. We stop. We need kjpindex*nbdl words = ',&
1251            & kjpindex*nbdl
1252       STOP 'sechiba_init'
1253    END IF
1254    ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
1255    IF (ier.NE.0) THEN
1256       WRITE (numout,*) ' error in co2_flux allocation. We stop. We need kjpindex*nvm words = ' ,kjpindex*nvm
1257       STOP 'sechiba_init'
1258    END IF
1259    co2_flux(:,:)=zero
1260
1261    ALLOCATE (shumdiag(kjpindex,nbdl),stat=ier)
1262    IF (ier.NE.0) THEN
1263       WRITE (numout,*) ' error in shumdiag allocation. We stop. We need kjpindex*nbdl words = ',&
1264            & kjpindex*nbdl
1265       STOP 'sechiba_init'
1266    END IF
1267
1268    ALLOCATE (litterhumdiag(kjpindex),stat=ier)
1269    IF (ier.NE.0) THEN
1270       WRITE (numout,*) ' error in litterhumdiag allocation. We stop. We need kjpindex words = ',kjpindex
1271       STOP 'sechiba_init'
1272    END IF
1273
1274    ! 1.4 two dimensions array allocation
1275
1276    ALLOCATE (vevapwet(kjpindex,nvm),stat=ier)
1277    IF (ier.NE.0) THEN
1278       WRITE (numout,*) ' error in vevapwet allocation. We stop. We need kjpindex x nvm words = ',&
1279            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1280       STOP 'sechiba_init'
1281    END IF
1282
1283    ALLOCATE (transpir(kjpindex,nvm),stat=ier)
1284    IF (ier.NE.0) THEN
1285       WRITE (numout,*) ' error in transpir allocation. We stop. We need kjpindex x nvm words = ',&
1286            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1287       STOP 'sechiba_init'
1288    END IF
1289
1290    ALLOCATE (qsintmax(kjpindex,nvm),stat=ier)
1291    IF (ier.NE.0) THEN
1292       WRITE (numout,*) ' error in qsintmax allocation. We stop. We need kjpindex x nvm words = ',&
1293            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm 
1294       STOP 'sechiba_init'
1295    END IF
1296
1297    ALLOCATE (rveget(kjpindex,nvm),stat=ier)
1298    IF (ier.NE.0) THEN
1299       WRITE (numout,*) ' error in rveget allocation. We stop. We need kjpindex x nvm words = ',&
1300            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1301       STOP 'sechiba_init'
1302    END IF
1303
1304    !
1305    !  1.5  Get the indexing table for the vegetation fields. In SECHIBA we work on reduced grids but to store in the
1306    !          full 3D filed vegetation variable we need another index table : indexveg, indexsoil, indexnobio and
1307    !          indexgrnd
1308    !
1309    DO ji = 1, kjpindex
1310       !
1311       DO jv = 1, nvm
1312          indexveg((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1313       ENDDO
1314       !     
1315       DO jv = 1, nstm
1316          indexsoil((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1317       ENDDO
1318       !     
1319       DO jv = 1, nnobio
1320          indexnobio((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1321       ENDDO
1322       !
1323       DO jv = 1, ngrnd
1324          indexgrnd((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1325       ENDDO
1326       !
1327       DO jv = 1, nslm
1328          indexlayer((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1329       ENDDO
1330       !
1331       DO jv = 1, 2
1332          indexalb((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
1333       ENDDO
1334       !
1335    ENDDO
1336
1337    !
1338    ! 2. restart value
1339    !
1340    ! open restart input file if needed
1341    ! and read data from restart input file
1342    !
1343
1344    IF (ldrestart_read) THEN
1345
1346       IF (long_print) WRITE (numout,*) ' we have to read a restart file for SECHIBA variables'
1347       !
1348       !   Read the default value that will be put into variable which are not in the restart file
1349       !
1350       CALL ioget_expval(val_exp)
1351
1352    ENDIF
1353
1354    !
1355    river_routing = control_in%river_routing
1356    hydrol_cwrr = control_in%hydrol_cwrr
1357    !
1358    ! 4. run control: store flags in a common variable
1359    !
1360
1361    control%river_routing = control_in%river_routing
1362    control%hydrol_cwrr = control_in%hydrol_cwrr
1363    control%ok_co2 = control_in%ok_co2
1364    control%ok_sechiba = control_in%ok_sechiba
1365    control%ok_stomate = control_in%ok_stomate
1366    control%ok_dgvm = control_in%ok_dgvm
1367    control%ok_pheno = control_in%ok_pheno
1368    control%stomate_watchout = control_in%stomate_watchout
1369
1370    IF (long_print) WRITE (numout,*) ' sechiba_init done '
1371
1372  END SUBROUTINE sechiba_init
1373  !
1374  !------------------------------------------------------------------
1375  !
1376  SUBROUTINE sechiba_clear (forcing_name,cforcing_name)
1377
1378    CHARACTER(LEN=100), INTENT(in)           :: forcing_name
1379    CHARACTER(LEN=100), INTENT(in)           :: cforcing_name
1380
1381    !
1382    ! initialisation
1383    !
1384    l_first_sechiba=.TRUE.
1385
1386    ! 1. Deallocate all dynamic variables
1387
1388    IF ( ALLOCATED (indexveg)) DEALLOCATE (indexveg)
1389    IF ( ALLOCATED (indexsoil)) DEALLOCATE (indexsoil)
1390    IF ( ALLOCATED (indexnobio)) DEALLOCATE (indexnobio)
1391    IF ( ALLOCATED (indexgrnd)) DEALLOCATE (indexgrnd)
1392    IF ( ALLOCATED (indexlayer)) DEALLOCATE (indexlayer)
1393    IF ( ALLOCATED (indexalb)) DEALLOCATE (indexalb)
1394    IF ( ALLOCATED (snow)) DEALLOCATE (snow)
1395    IF ( ALLOCATED (snow_age)) DEALLOCATE (snow_age)
1396    IF ( ALLOCATED (drysoil_frac)) DEALLOCATE (drysoil_frac)
1397    IF ( ALLOCATED (rsol)) DEALLOCATE (rsol)
1398    IF ( ALLOCATED (evap_bare_lim)) DEALLOCATE (evap_bare_lim)
1399    IF ( ALLOCATED (evapot)) DEALLOCATE (evapot)
1400    IF ( ALLOCATED (evapot_corr)) DEALLOCATE (evapot_corr)
1401    IF ( ALLOCATED (humrel)) DEALLOCATE (humrel)
1402    IF ( ALLOCATED (vegstress)) DEALLOCATE (vegstress)
1403    IF ( ALLOCATED (soiltype)) DEALLOCATE (soiltype)
1404    IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1)
1405    IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4)
1406    IF ( ALLOCATED (soilcap)) DEALLOCATE (soilcap)
1407    IF ( ALLOCATED (soilflx)) DEALLOCATE (soilflx)
1408    IF ( ALLOCATED (temp_sol)) DEALLOCATE (temp_sol)
1409    IF ( ALLOCATED (qsurf)) DEALLOCATE (qsurf)
1410    IF ( ALLOCATED (qsintveg)) DEALLOCATE (qsintveg)
1411    IF ( ALLOCATED (vbeta2))  DEALLOCATE (vbeta2)
1412    IF ( ALLOCATED (vbeta3)) DEALLOCATE (vbeta3)
1413    IF ( ALLOCATED (vbetaco2)) DEALLOCATE (vbetaco2)
1414    IF ( ALLOCATED (cimean)) DEALLOCATE (cimean)
1415    IF ( ALLOCATED (gpp)) DEALLOCATE (gpp)
1416    IF ( ALLOCATED (veget)) DEALLOCATE (veget)
1417    IF ( ALLOCATED (veget_max)) DEALLOCATE (veget_max)
1418    IF ( ALLOCATED (lai)) DEALLOCATE (lai)
1419    IF ( ALLOCATED (height)) DEALLOCATE (height)
1420    IF ( ALLOCATED (roughheight)) DEALLOCATE (roughheight)
1421    IF ( ALLOCATED (frac_nobio)) DEALLOCATE (frac_nobio)
1422    IF ( ALLOCATED (snow_nobio)) DEALLOCATE (snow_nobio)
1423    IF ( ALLOCATED (snow_nobio_age)) DEALLOCATE (snow_nobio_age)
1424    IF ( ALLOCATED (assim_param)) DEALLOCATE (assim_param)
1425    IF ( ALLOCATED (vevapsno)) DEALLOCATE (vevapsno)
1426    IF ( ALLOCATED (vevapnu)) DEALLOCATE (vevapnu)
1427    IF ( ALLOCATED (t2mdiag)) DEALLOCATE (t2mdiag)
1428    IF ( ALLOCATED (totfrac_nobio)) DEALLOCATE (totfrac_nobio)
1429    IF ( ALLOCATED (runoff)) DEALLOCATE (runoff)
1430    IF ( ALLOCATED (drainage)) DEALLOCATE (drainage)
1431    IF ( ALLOCATED (returnflow)) DEALLOCATE (returnflow)
1432    IF ( ALLOCATED (irrigation)) DEALLOCATE (irrigation)
1433    IF ( ALLOCATED (tot_melt)) DEALLOCATE (tot_melt)
1434    IF ( ALLOCATED (valpha)) DEALLOCATE (valpha)
1435    IF ( ALLOCATED (vbeta)) DEALLOCATE (vbeta)
1436    IF ( ALLOCATED (fusion)) DEALLOCATE (fusion)
1437    IF ( ALLOCATED (rau)) DEALLOCATE (rau)
1438    IF ( ALLOCATED (deadleaf_cover)) DEALLOCATE (deadleaf_cover)
1439    IF ( ALLOCATED (stempdiag)) DEALLOCATE (stempdiag)
1440    IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
1441    IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
1442    IF ( ALLOCATED (litterhumdiag)) DEALLOCATE (litterhumdiag)
1443    IF ( ALLOCATED (vevapwet)) DEALLOCATE (vevapwet)
1444    IF ( ALLOCATED (transpir)) DEALLOCATE (transpir)
1445    IF ( ALLOCATED (qsintmax)) DEALLOCATE (qsintmax)
1446    IF ( ALLOCATED (rveget)) DEALLOCATE (rveget)
1447    ! 2. clear all modules
1448    CALL pft_parameters_clear
1449    CALL slowproc_clear 
1450    CALL diffuco_clear 
1451    CALL enerbil_clear 
1452    IF ( hydrol_cwrr ) THEN
1453       CALL hydrol_clear 
1454    ELSE
1455       CALL hydrolc_clear 
1456    ENDIF
1457    CALL condveg_clear 
1458    CALL thermosoil_clear
1459    CALL routing_clear
1460    !3. give name to next block
1461    stomate_forcing_name=forcing_name
1462    stomate_Cforcing_name=Cforcing_name
1463
1464  END SUBROUTINE sechiba_clear
1465
1466  !! SECHIBA's variables initialisation
1467  !! called every time step
1468  !!
1469  SUBROUTINE sechiba_var_init (kjpindex, rau, pb, temp_air) 
1470
1471    ! interface description
1472    ! input scalar
1473    INTEGER(i_std), INTENT (in)                    :: kjpindex           !! Domain dimension
1474    ! input fields
1475    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: pb                 !! Lowest level pressure
1476    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: temp_air           !! Air temperature
1477    ! output fields
1478    REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rau                !! Density
1479
1480    ! local declaration
1481    INTEGER(i_std)                                :: ji
1482
1483    !
1484    ! initialisation
1485    !
1486
1487    !
1488    ! 1. calcul of rau: air density
1489    !
1490
1491    DO ji = 1,kjpindex
1492       rau(ji) = pa_par_hpa * pb(ji) / (cte_molr*temp_air(ji))
1493    END DO
1494
1495    IF (long_print) WRITE (numout,*) ' sechiba_var_init done '
1496
1497  END SUBROUTINE sechiba_var_init
1498
1499  !!
1500  !! Swap new fields to previous fields
1501  !!
1502  SUBROUTINE sechiba_end (kjpindex, dtradia, temp_sol, temp_sol_new)
1503
1504    ! interface description
1505    ! input scalar
1506    INTEGER(i_std), INTENT (in)                       :: kjpindex           !! Domain dimension
1507    REAL(r_std),INTENT (in)                           :: dtradia            !! Time step in seconds
1508    ! input fields
1509    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: temp_sol_new       !! New soil temperature
1510    ! output fields
1511    REAL(r_std),DIMENSION (kjpindex), INTENT (out)    :: temp_sol           !! Soil temperature
1512
1513    !
1514    ! swap
1515    !
1516    temp_sol(:) = temp_sol_new(:)
1517
1518    IF (long_print) WRITE (numout,*) ' sechiba_end done '
1519
1520  END SUBROUTINE sechiba_end
1521
1522END MODULE sechiba
Note: See TracBrowser for help on using the repository browser.