source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_sechiba/sechiba.f90 @ 6268

Last change on this file since 6268 was 734, checked in by didier.solyga, 12 years ago

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