source: branches/publications/ORCHIDEE-Clateral/src_sechiba/sechiba.f90 @ 7329

Last change on this file since 7329 was 7191, checked in by haicheng.zhang, 3 years ago

Implementing latral transfers of sediment and POC from land to ocean through river in ORCHILEAK

  • Property svn:keywords set to Revision Date HeadURL Date Author Revision
File size: 155.4 KB
Line 
1!  ==============================================================================================================================\n
2!  MODULE       : sechiba
3!
4!  CONTACT      : orchidee-help _at_ listes.ipsl.fr
5!
6!  LICENCE      : IPSL (2006)
7!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8!
9!>\BRIEF        Structures the calculation of atmospheric and hydrological
10!! variables by calling diffuco_main, enerbil_main, hydrolc_main (or hydrol_main),
11!! enerbil_fusion, condveg_main and thermosoil_main. Note that sechiba_main
12!! calls slowproc_main and thus indirectly calculates the biogeochemical
13!! processes as well.
14!!
15!!\n DESCRIPTION  : :: shumdiag, :: litterhumdiag and :: stempdiag are not
16!! saved in the restart file because at the first time step because they
17!! are recalculated. However, they must be saved as they are in slowproc
18!! which is called before the modules which calculate them.
19!!
20!! RECENT CHANGE(S): None
21!!
22!! REFERENCE(S) : None
23!!   
24!! SVN     :
25!! $HeadURL$
26!! $Date$
27!! $Revision$
28!! \n
29!_ ================================================================================================================================
30 
31MODULE sechiba
32 
33  USE ioipsl
34  USE xios_orchidee
35  USE constantes
36  USE time, ONLY : one_day, dt_sechiba
37  USE constantes_soil
38  USE pft_parameters
39  USE grid
40  USE diffuco
41  USE condveg
42  USE enerbil
43  USE hydrol
44  USE hydrolc
45  USE thermosoil
46  USE thermosoilc
47  USE sechiba_io_p
48  USE slowproc
49  USE erosion
50  USE routing
51  USE ioipsl_para
52  USE chemistry
53
54  IMPLICIT NONE
55
56  PRIVATE
57  PUBLIC sechiba_main, sechiba_initialize, sechiba_clear, sechiba_interface_orchidee_inca
58
59  INTEGER(i_std), SAVE                             :: printlev_loc   !! local printlev for this module
60!$OMP THREADPRIVATE(printlev_loc)
61  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexveg       !! indexing array for the 3D fields of vegetation
62!$OMP THREADPRIVATE(indexveg)
63  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlai       !! indexing array for the 3D fields of vegetation
64!$OMP THREADPRIVATE(indexlai)
65  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnobio     !! indexing array for the 3D fields of other surfaces (ice,
66                                                                     !! lakes, ...)
67!$OMP THREADPRIVATE(indexnobio)
68  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsoil      !! indexing array for the 3D fields of soil types (kjpindex*nstm)
69!$OMP THREADPRIVATE(indexsoil)
70  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexgrnd      !! indexing array for the 3D ground heat profiles (kjpindex*ngrnd)
71!$OMP THREADPRIVATE(indexgrnd)
72  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlayer     !! indexing array for the 3D fields of soil layers in CWRR (kjpindex*nslm)
73!$OMP THREADPRIVATE(indexlayer)
74  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnslm      !! indexing array for the 3D fields of diagnostic soil layers (kjpindex*nslm)
75!$OMP THREADPRIVATE(indexnslm)
76  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexalb       !! indexing array for the 2 fields of albedo
77!$OMP THREADPRIVATE(indexalb)
78  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsnow      !! indexing array for the 3D fields snow layers
79!$OMP THREADPRIVATE(indexsnow)
80  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: veget          !! Fraction of vegetation type (unitless, 0-1)       
81!$OMP THREADPRIVATE(veget)
82  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: veget_max      !! Max. fraction of vegetation type (LAI -> infty, unitless)
83!$OMP THREADPRIVATE(veget_max)
84  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: tot_bare_soil  !! Total evaporating bare soil fraction
85!$OMP THREADPRIVATE(tot_bare_soil)
86  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: height         !! Vegetation Height (m)
87!$OMP THREADPRIVATE(height)
88  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: totfrac_nobio  !! Total fraction of continental ice+lakes+cities+...
89                                                                     !! (unitless, 0-1)
90!$OMP THREADPRIVATE(totfrac_nobio)
91  REAL(r_std),ALLOCATABLE, SAVE,  DIMENSION (:,:)  :: wat_flux0      !! Water flux in the first soil layers exported for soil C calculations
92!$OMP THREADPRIVATE(wat_flux0)
93  REAL(r_std),ALLOCATABLE, SAVE,  DIMENSION (:,:,:) :: wat_flux       !! Water flux in the soil layers exported for soil C calculations
94!$OMP THREADPRIVATE(wat_flux)
95  REAL(r_std),ALLOCATABLE, SAVE,  DIMENSION (:,:)  :: drainage_per_soil !! Drainage per soil type exported for soil C calculations
96!$OMP THREADPRIVATE(drainage_per_soil)
97  REAL(r_std),ALLOCATABLE, SAVE,  DIMENSION (:,:)  :: runoff_per_soil !! Runoff per soil type exported for soil C calculations
98!$OMP THREADPRIVATE(runoff_per_soil)
99  REAL, ALLOCATABLE, DIMENSION(:,:,:) :: soil_mc
100!$OMP THREADPRIVATE(soil_mc)
101  REAL, ALLOCATABLE, DIMENSION(:,:) :: litter_mc
102!$OMP THREADPRIVATE(litter_mc)                                                                                                                                                                                                                                                                                 
103 !!
104!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Added by Haicheng Zhang(ZHC) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
105  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: laydep           !! The lower limit of each (11) soil layer (m) ::zhc
106  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: bulkdens         !! bulk_density (kg m-3)    zhc
107  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)       :: textfrac         !! Fraction of different soil particle sizes (0-1)
108  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: gridarea         !! grid area (m^2)
109  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: effgrid_perc     !! Percentage of effective area in an ORCHIDEE pixel (0-1)
110  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: std_totsed       !! Standrad sediment loss rate(ton)    zhc
111  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: totrunoff_d      !! daily total runoff (mm) 
112  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: peakrunoff_d     !! daily peak runoff rate.(mm/s = kg/m-2/s)
113  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: totrundra_d      !! daily total runoff (mm) 
114  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: peakrundra_d     !! daily peak runoff rate.(mm/s = kg/m-2/s)
115  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)       :: erocap           !! Erosion capacity(t/day)                                                                                                                           
116  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)       :: erodep           !! Depth of the eroded soil under different PFTs (m)
117  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)     :: eroc                 !! Eroded SOC(gC/m^2/day) for each PFT
118  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)               :: POC_EXP_agg          !! POC exports along with surface runoff, caused by soil erosion, in 
119                                                                                                                                                        !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
120  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)               :: DOC_ero_agg          !! DOC exports along with surface runoff,caused by soil erosion, in 
121                                                                                                                                                        !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
122  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)               :: SED_EXP_agg          !! sediment exports along with surface runoff,caused by soil erosion, in 
123                                                                                                                                                        !! @tex $(g m^{-2} dt_sechiba^{-1})$ @endtex
124  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: tot_erocap       !! Erosion capacity(t/day)                                                                                                                           
125  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)         :: ave_erodep       !! Depth of the eroded soil in orchidee grids (m)
126  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)       :: tot_eroc         !! Eroded SOC(gC/m**2)
127!  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)     :: eron                !! Eroded SOC(gC/m**2)
128  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)        :: cf_veget         !! Daily vegetation cover factor
129  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)        :: kf_litter        !! Daily above-ground litter cover factor
130  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)        :: kf_root          !! Daily above-ground litter cover factor
131  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:)       :: rootmass         !! belowground biomass
132                                                                            !! per ground area                                                                             !! @tex $(gC m^{-2})$ @endtex
133  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:)       :: litter_above     !! Above ground metabolic and structural litter
134                                                                            !! per ground area
135                                                                            !! @tex $(gC m^{-2})$ @endtex
136  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:)     :: litter_below     !! Below ground metabolic and structural litter
137                                                                            !! per ground area
138                                                                            !! @tex $(gC m^{-2})$ @endtex
139  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:)       :: carbon           !! Soil carbon pools per ground area: active, slow, or
140                                                                            !! passive, @tex $(gC m^{-2})$ @endtex
141  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:,:)   :: DOC              !! Soil dissolved organic carbon free or adsorbed
142                                                                            !! detailled for each pools @tex $(gC m^{-2} of ground)$ @endtex
143  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)           :: lignin_struc_above   !! Ratio Lignine/Carbon in structural litter for above
144                                                                            !! ground compartments (unitless)
145  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)         :: lignin_struc_below   !! Ratio Lignine/Carbon in structural litter for below
146                                                                            !! ground compartments (unitless)                                                                                                                                   
147  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)        :: depth_deepsoil   !! Depth of the soil layer deeper than 2 m.
148                                                                            !! When sediment deposition occuring, the original surface (0-2)
149                                                                                                                                                        !! soil DOC, and SOC will enther into this layer.
150  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)        :: seddep_rate      !! Sediment deposition rate (m day-1) for each PFT
151                                                                            !! When sediment deposition occuring, the original surface (0-2)
152                                                                                                                                                        !! soil DOC, and SOC will enther into this layer.                                                                                                                                                       
153  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)          :: sed_deposition   !! Mean sediment deposition rate at
154                                                                            !! floodplains (for each sechiba time step, dt_sechiba)
155                                                                            !! @tex $(kg m^{-2} 30-munites-1)$ @endtex
156                                                                                                                                                       
157  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)        :: poc_deposition   !! Mean POC deposition rate at
158                                                                            !! floodplains (for each sechiba time step, dt_sechiba)
159                                                                            !! @tex $(kg C m^{-2} 30-munites-1)$ @endtex
160  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)        :: sed_deposition_d !! Mean daily sediment deposition rate at
161                                                                            !! floodplains (for each sechiba time step, dt_sechiba)
162                                                                            !! @tex $(kg m^{-2} day-1)$ @endtex
163                                                                                                                                                       
164  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)        :: poc_deposition_d !! Mean daily POC deposition rate at
165                                                                            !! floodplains (for each sechiba time step, dt_sechiba)
166                                                                            !! @tex $(kg C m^{-2} day-1)$ @endtex                                                                                                                                                       
167!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Added by Haicheng Zhang(ZHC) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
168!!
169  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: floodout       !! Flow out of floodplains from hydrol
170!$OMP THREADPRIVATE(floodout)
171  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: runoff         !! Surface runoff calculated by hydrol or hydrolc
172                                                                     !! @tex $(kg m^{-2})$ @endtex
173  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: DOC_EXP_agg    !! DOC exports, diffrenet paths (nexp), in 
174                                                                     !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
175!$OMP THREADPRIVATE(runoff)
176  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: drainage       !! Deep drainage calculatedd by hydrol or hydrolc
177                                                                     !! @tex $(kg m^{-2})$ @endtex
178                                                                     !! @tex $(kg m^{-2})$ @endtex
179!$OMP THREADPRIVATE(drainage)
180  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: returnflow     !! Water flow from lakes and swamps which returns to
181                                                                     !! the grid box @tex $(kg m^{-2})$ @endtex
182!$OMP THREADPRIVATE(returnflow)
183  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: reinfiltration !! Routed water which returns into the soil
184!$OMP THREADPRIVATE(reinfiltration)
185  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: irrigation     !! Irrigation flux taken from the routing reservoirs and
186                                                                     !! being put into the upper layers of the soil
187                                                                     !! @tex $(kg m^{-2})$ @endtex                                                                                                                               
188!$OMP THREADPRIVATE(irrigation)
189  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: DOC_to_topsoil !! Sum of DOC fluxes from irrirgation and reinfiltration
190                                                                     !! @tex $(g m^{-2})$ @endtex                                           
191!$OMP THREADPRIVATE(DOC_to_topsoil)
192  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: DOC_to_subsoil !! DOC fluxes from returnflow in lakes and swamps
193                                                                     !! @tex $(g m^{-2})$ @endtex                                           
194!$OMP THREADPRIVATE(DOC_to_subsoil)
195  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)      :: precip2canopy  !! Precipitation onto the canopy
196!$OMP THREADPRIVATE(precip2canopy)
197  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)      :: precip2ground  !! Precipitation not intercepted by canopy
198!$OMP THREADPRIVATE(precip2ground)
199  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)      :: canopy2ground  !! Water flux from canopy to the ground
200!$OMP THREADPRIVATE(canopy2ground) 
201  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: emis           !! Surface emissivity (unitless)
202!$OMP THREADPRIVATE(emis)
203  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: z0h           !! Surface roughness for heat (m)
204!$OMP THREADPRIVATE(z0h)
205  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: z0m           !! Surface roughness for momentum (m)
206!$OMP THREADPRIVATE(z0m)
207  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: roughheight    !! Effective height for roughness (m)
208!$OMP THREADPRIVATE(roughheight)
209  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: reinf_slope    !! slope coefficient (reinfiltration)
210!$OMP THREADPRIVATE(reinf_slope)
211  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: shumdiag       !! Mean relative soil moisture in the different levels used
212                                                                     !! by thermosoil.f90 (unitless, 0-1)
213!$OMP THREADPRIVATE(shumdiag)
214  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: shumdiag_perma !! Saturation degree of the soil
215!$OMP THREADPRIVATE(shumdiag_perma)
216  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: k_litt         !! litter cond.
217!$OMP THREADPRIVATE(k_litt)
218  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: litterhumdiag  !! Litter dryness factor (unitless, 0-1)
219!$OMP THREADPRIVATE(litterhumdiag)
220  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: stempdiag      !! Temperature which controls canopy evolution (K)
221!$OMP THREADPRIVATE(stempdiag)
222  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: qsintveg       !! Water on vegetation due to interception
223                                                                     !! @tex $(kg m^{-2})$ @endtex
224!$OMP THREADPRIVATE(qsintveg)
225  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta2         !! Interception resistance (unitless,0-1)
226!$OMP THREADPRIVATE(vbeta2)
227  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta3         !! Vegetation resistance (unitless,0-1)
228!$OMP THREADPRIVATE(vbeta3)
229  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta3pot      !! Potential vegetation resistance
230!$OMP THREADPRIVATE(vbeta3pot)
231  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: gsmean         !! Mean stomatal conductance for CO2 (mol m-2 s-1)
232!$OMP THREADPRIVATE(gsmean)
233  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: cimean         !! STOMATE: mean intercellular CO2 concentration (ppm)
234!$OMP THREADPRIVATE(cimean)
235  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vevapwet       !! Interception loss over each PFT
236                                                                     !! @tex $(kg m^{-2} days^{-1})$ @endtex
237!$OMP THREADPRIVATE(vevapwet)
238  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: transpir       !! Transpiration @tex $(kg m^{-2} days^{-1})$ @endtex
239!$OMP THREADPRIVATE(transpir)
240  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: transpot       !! Potential Transpiration (needed for irrigation)
241!$OMP THREADPRIVATE(transpot)
242  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: qsintmax       !! Maximum amount of water in the canopy interception
243                                                                     !! reservoir @tex $(kg m^{-2})$ @endtex
244!$OMP THREADPRIVATE(qsintmax)
245  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: rveget         !! Surface resistance for the vegetation
246                                                                     !! @tex $(s m^{-1})$ @endtex
247!$OMP THREADPRIVATE(rveget)
248  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: rstruct        !! Vegetation structural resistance
249!$OMP THREADPRIVATE(rstruct)
250  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: snow_nobio     !! Snow mass of non-vegetative surfaces
251                                                                     !! @tex $(kg m^{-2})$ @endtex
252!$OMP THREADPRIVATE(snow_nobio)
253  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: snow_nobio_age !! Snow age on non-vegetative surfaces (days)
254!$OMP THREADPRIVATE(snow_nobio_age)
255  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: frac_nobio     !! Fraction of non-vegetative surfaces (continental ice,
256                                                                     !! lakes, ...) (unitless, 0-1)
257!$OMP THREADPRIVATE(frac_nobio)
258  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: assim_param    !! min+max+opt temps, vcmax, vjmax for photosynthesis
259!$OMP THREADPRIVATE(assim_param)
260  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lai            !! Surface foliaire
261!$OMP THREADPRIVATE(lai)
262  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: gpp            !! STOMATE: GPP. gC/m**2 of total area
263!$OMP THREADPRIVATE(gpp)
264  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)     :: temp_growth    !! Growth temperature (C) - Is equal to t2m_month
265!$OMP THREADPRIVATE(temp_growth)
266  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: humrel         !! Relative humidity
267!$OMP THREADPRIVATE(humrel)
268  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vegstress      !! Vegetation moisture stress (only for vegetation growth)
269!$OMP THREADPRIVATE(vegstress)
270  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: frac_age       !! Age efficacity from STOMATE for isoprene
271!$OMP THREADPRIVATE(frac_age)
272  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: soiltile       !! Fraction of each soil tile (0-1, unitless)
273!$OMP THREADPRIVATE(soiltile)
274  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: fraclut        !! Fraction of each landuse tile (0-1, unitless)
275!$OMP THREADPRIVATE(fraclut)
276  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: nwdFraclut     !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless)
277!$OMP THREADPRIVATE(nwdFraclut)
278  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: njsc           !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
279!$OMP THREADPRIVATE(njsc)
280  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta1         !! Snow resistance
281!$OMP THREADPRIVATE(vbeta1)
282  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta4         !! Bare soil resistance
283!$OMP THREADPRIVATE(vbeta4)
284  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta5         !! Floodplains resistance
285!$OMP THREADPRIVATE(vbeta5)
286  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: soilcap        !!
287!$OMP THREADPRIVATE(soilcap)
288  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: soilflx        !!
289!$OMP THREADPRIVATE(soilflx)
290  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: temp_sol       !! Soil temperature
291!$OMP THREADPRIVATE(temp_sol)
292  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: qsurf          !! near soil air moisture
293!$OMP THREADPRIVATE(qsurf)
294  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: flood_res      !! flood reservoir estimate
295!$OMP THREADPRIVATE(flood_res)
296  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: fastr          !! fast reservoir estimate
297!$OMP THREADPRIVATE(fastr)
298  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: flood_frac     !! flooded fraction
299!$OMP THREADPRIVATE(flood_frac)
300  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: stream_frac     !! stream fraction
301!$OMP THREADPRIVATE(stream_frac)
302  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: streamfl_frac   !! flooded by swollen stream fraction
303!$OMP THREADPRIVATE(streamfl_frac)
304  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: snow           !! Snow mass [Kg/m^2]
305!$OMP THREADPRIVATE(snow)
306  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: snow_age       !! Snow age @tex ($d$) @endtex
307!$OMP THREADPRIVATE(snow_age)
308  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: drysoil_frac   !! Fraction of visibly (albedo) Dry soil (Between 0 and 1)
309!$OMP THREADPRIVATE(drysoil_frac)
310  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: rsol           !! resistance to bare soil evaporation
311!$OMP THREADPRIVATE(rsol)
312  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evap_bare_lim  !! Bare soil stress
313!$OMP THREADPRIVATE(evap_bare_lim)
314
315  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: co2_flux       !! CO2 flux (gC/m**2 of average ground/s)
316!$OMP THREADPRIVATE(co2_flux)
317  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: co2_to_bm      !! virtual CO2 flux (gC/m**2 of average ground/s)
318!$OMP THREADPRIVATE(co2_to_bm)
319  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evapot         !! Soil Potential Evaporation
320!$OMP THREADPRIVATE(evapot)
321  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evapot_corr    !! Soil Potential Evaporation Correction (Milly 1992)
322!$OMP THREADPRIVATE(evapot_corr)
323  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapflo       !! Floodplains evaporation
324!$OMP THREADPRIVATE(vevapflo)
325  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapsno       !! Snow evaporation
326!$OMP THREADPRIVATE(vevapsno)
327  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapnu        !! Bare soil evaporation
328!$OMP THREADPRIVATE(vevapnu)
329  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: tot_melt       !! Total melt
330!$OMP THREADPRIVATE(tot_melt)
331  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta          !! Resistance coefficient
332!$OMP THREADPRIVATE(vbeta)
333  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: fusion         !!
334!$OMP THREADPRIVATE(fusion)
335  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: rau            !! Density
336!$OMP THREADPRIVATE(rau)
337  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: deadleaf_cover !! Fraction of soil covered by dead leaves
338!$OMP THREADPRIVATE(deadleaf_cover)
339  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: ptnlev1        !! 1st level Different levels soil temperature
340!$OMP THREADPRIVATE(ptnlev1)
341  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: mc_layh        !! Volumetric soil moisture for each layer in hydrol(liquid + ice) (m3/m3)
342!$OMP THREADPRIVATE(mc_layh)
343  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: mcl_layh       !! Volumetric soil moisture for each layer in hydrol(liquid) (m3/m3)
344!$OMP THREADPRIVATE(mcl_layh)
345  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: soilmoist      !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
346!$OMP THREADPRIVATE(soilmoist)
347
348  LOGICAL, SAVE                                    :: l_first_sechiba = .TRUE. !! Flag controlling the intialisation (true/false)
349!$OMP THREADPRIVATE(l_first_sechiba)
350
351  ! Variables related to snow processes calculations 
352
353  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: frac_snow_veg   !! Snow cover fraction on vegetation (unitless)
354!$OMP THREADPRIVATE(frac_snow_veg)
355  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: frac_snow_nobio !! Snow cover fraction on continental ice, lakes, etc (unitless)
356!$OMP THREADPRIVATE(frac_snow_nobio)
357  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowrho      !! snow density for each layer
358!$OMP THREADPRIVATE(snowrho)
359  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowheat     !! snow heat content for each layer (J/m2)
360!$OMP THREADPRIVATE(snowheat)
361  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowgrain    !! snow grain size (m)
362!$OMP THREADPRIVATE(snowgrain)
363  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowtemp     !! snow temperature profile (K)
364!$OMP THREADPRIVATE(snowtemp)
365  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowdz       !! snow layer thickness (m)
366!$OMP THREADPRIVATE(snowdz)
367  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: gtemp        !! soil surface temperature
368!$OMP THREADPRIVATE(gtemp)
369  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: pgflux       !! net energy into snow pack
370!$OMP THREADPRIVATE(pgflux)
371  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: cgrnd_snow   !! Integration coefficient for snow numerical scheme
372!$OMP THREADPRIVATE(cgrnd_snow)
373  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: dgrnd_snow   !! Integration coefficient for snow numerical scheme
374!$OMP THREADPRIVATE(dgrnd_snow)
375  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: lambda_snow  !! Coefficient of the linear extrapolation of surface temperature
376                                                                  !! from the first and second snow layers
377!$OMP THREADPRIVATE(lambda_snow)
378  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
379!$OMP THREADPRIVATE(temp_sol_add)
380
381  LOGICAL, SAVE                                   :: ok_doc= .TRUE.
382!$OMP THREADPRIVATE(ok_doc)   
383  CHARACTER(LEN=80)    :: var_name       !! To store variables names for I/O (unitless)
384  CHARACTER(LEN=10)    :: flow_suff      !! To assign a suffix to the variables name indicating the transported specie
385  CHARACTER(LEN=80)    :: var_long_name  !! To store variables long names for I/O (unitless)
386CONTAINS
387
388!!  =============================================================================================================================
389!! SUBROUTINE:    sechiba_initialize
390!!
391!>\BRIEF          Initialize all prinicipal modules by calling their "_initialize" subroutines
392!!
393!! DESCRIPTION:   Initialize all prinicipal modules by calling their "_initialize" subroutines
394!!
395!! \n
396!_ ==============================================================================================================================
397
398  SUBROUTINE sechiba_initialize( &
399       kjit,         kjpij,        kjpindex,     index,                   &
400       lalo,         contfrac,     neighbours,   resolution, zlev,        &
401       u,            v,            qair,         t2m,        temp_air,    &
402       petAcoef,     peqAcoef,     petBcoef,     peqBcoef,                &
403       precip_rain,  precip_snow,  lwdown,       swnet,      swdown,      &
404       pb,           rest_id,      hist_id,      hist2_id,                &
405       rest_id_stom, hist_id_stom, hist_id_stom_IPCC,                     &
406       coastalflow,  riverflow,    tsol_rad,     vevapp,       qsurf_out, &
407       z0m_out,      z0h_out,      albedo,       fluxsens,     fluxlat,      emis_out,  &
408       netco2flux,   fco2_lu,      temp_sol_new, tq_cdrag)
409
410!! 0.1 Input variables
411    INTEGER(i_std), INTENT(in)                               :: kjit              !! Time step number (unitless)
412    INTEGER(i_std), INTENT(in)                               :: kjpij             !! Total size of the un-compressed grid
413                                                                                  !! (unitless)
414    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size - terrestrial pixels only
415                                                                                  !! (unitless)
416    INTEGER(i_std),INTENT (in)                               :: rest_id           !! _Restart_ file identifier (unitless)
417    INTEGER(i_std),INTENT (in)                               :: hist_id           !! _History_ file identifier (unitless)
418    INTEGER(i_std),INTENT (in)                               :: hist2_id          !! _History_ file 2 identifier (unitless)
419    INTEGER(i_std),INTENT (in)                               :: rest_id_stom      !! STOMATE's _Restart_ file identifier
420                                                                                  !! (unitless)
421    INTEGER(i_std),INTENT (in)                               :: hist_id_stom      !! STOMATE's _History_ file identifier
422                                                                                  !! (unitless)
423    INTEGER(i_std),INTENT(in)                                :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file
424                                                                                  !! identifier (unitless)
425    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)          :: lalo              !! Geographic coordinates (latitude,longitude)
426                                                                                  !! for grid cells (degrees)
427    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: contfrac          !! Fraction of continent in the grid
428                                                                                  !! (unitless, 0-1)
429    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index             !! Indices of the pixels on the map.
430                                                                                  !! Sechiba uses a reduced grid excluding oceans
431                                                                                  !! ::index contains the indices of the
432                                                                                  !! terrestrial pixels only! (unitless)
433    INTEGER(i_std), DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours        !! Neighboring grid points if land!(unitless)
434    REAL(r_std), DIMENSION (kjpindex,2), INTENT(in)          :: resolution        !! Size in x and y of the grid (m)
435   
436    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u                 !! Lowest level wind speed in direction u
437                                                                                  !! @tex $(m.s^{-1})$ @endtex
438    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: v                 !! Lowest level wind speed in direction v
439                                                                                  !! @tex $(m.s^{-1})$ @endtex
440    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: zlev              !! Height of first layer (m)
441    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair              !! Lowest level specific humidity
442                                                                                  !! @tex $(kg kg^{-1})$ @endtex
443    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: t2m               !! 2m air temperature (K)
444    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: precip_rain       !! Rain precipitation
445                                                                                  !! @tex $(kg m^{-2})$ @endtex
446    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: precip_snow       !! Snow precipitation
447                                                                                  !! @tex $(kg m^{-2})$ @endtex
448    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: lwdown            !! Down-welling long-wave flux
449                                                                                  !! @tex $(W m^{-2})$ @endtex
450    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swnet             !! Net surface short-wave flux
451                                                                                  !! @tex $(W m^{-2})$ @endtex
452    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swdown            !! Down-welling surface short-wave flux
453                                                                                  !! @tex $(W m^{-2})$ @endtex
454    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_air          !! Air temperature (K)
455    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petAcoef          !! Coefficients A for T from the Planetary
456                                                                                  !! Boundary Layer
457    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqAcoef          !! Coefficients A for q from the Planetary
458                                                                                  !! Boundary Layer
459    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petBcoef          !! Coefficients B for T from the Planetary
460                                                                                  !! Boundary Layer
461    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqBcoef          !! Coefficients B for q from the Planetary
462                                                                                  !! Boundary Layer
463    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: pb                !! Surface pressure (hPa)
464
465
466!! 0.2 Output variables
467    REAL(r_std),DIMENSION (kjpindex,nflow), INTENT (out)     :: coastalflow       !! Outflow on coastal points by small basins.
468                                                                                  !! This is the water which flows in a disperse
469                                                                                  !! way into the ocean
470                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
471    REAL(r_std),DIMENSION (kjpindex,nflow), INTENT (out)     :: riverflow         !! Outflow of the major rivers.
472                                                                                  !! The flux will be located on the continental
473                                                                                  !! grid but this should be a coastal point 
474                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
475    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: tsol_rad          !! Radiative surface temperature
476                                                                                  !! @tex $(W m^{-2})$ @endtex
477    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapp            !! Total of evaporation
478                                                                                  !! @tex $(kg m^{-2} days^{-1})$ @endtex
479   
480    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qsurf_out         !! Surface specific humidity
481                                                                                  !! @tex $(kg kg^{-1})$ @endtex
482    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0m_out           !! Surface roughness momentum (output diagnostic, m)
483    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0h_out           !! Surface roughness heat (output diagnostic, m)
484    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out)         :: albedo            !! Surface albedo for visible and near-infrared
485                                                                                  !! (unitless, 0-1)
486    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxsens          !! Sensible heat flux
487                                                                                  !! @tex $(W m^{-2})$ @endtex
488    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxlat           !! Latent heat flux
489                                                                                  !! @tex $(W m^{-2})$ @endtex
490    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: emis_out          !! Emissivity (unitless)
491    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: netco2flux        !! Sum CO2 flux over PFTs
492                                                                                  !! ??(gC m^{-2} s^{-1})??
493    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fco2_lu           !! Land Cover Change CO2 flux
494                                                                                  !! ??(gC m^{-2} s^{-1})??
495    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: temp_sol_new      !! New ground temperature (K)
496
497!! 0.3 Modified
498    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: tq_cdrag          !! Surface drag coefficient (-)
499
500!! 0.4 Local variables
501    INTEGER(i_std)                                           :: ji, jv            !! Index (unitless)
502    REAL(r_std), DIMENSION(kjpindex)                         :: zmaxh_glo         !! 2D field of constant soil depth (zmaxh) (m)
503    CHARACTER(LEN=80)                                        :: var_name          !! To store variables names for I/O (unitless)
504
505!_ ================================================================================================================================
506
507    IF (printlev>=3) WRITE(numout,*) 'Start sechiba_initialize'
508
509    !! 1. Initialize variables on first call
510
511    !! 1.1 Initialize most of sechiba's variables
512    CALL sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)
513   
514    !! 1.3 Initialize stomate's variables
515        !WRITE(numout,*) 'Sechiba_ini_ZHC1'
516        !WRITE(numout,*) 'litter_above: ', SUM(litter_above(:,:,:,icarbon))
517        !WRITE(numout,*) 'litter_below: ', SUM(litter_below(:,:,:,:,icarbon))
518        !WRITE(numout,*) 'lignin_be',SUM(lignin_struc_below)
519        !WRITE(numout,*) 'DOC_SOC',SUM(DOC),SUM(carbon)
520   
521    CALL slowproc_initialize (kjit,          kjpij,        kjpindex,                          &
522                              rest_id,       rest_id_stom, hist_id_stom,   hist_id_stom_IPCC, &
523                              index,         indexveg,     lalo,           neighbours,        &
524                              resolution,    contfrac,     t2m,                               &
525                              soiltile,      reinf_slope,  deadleaf_cover, assim_param,       &
526                              lai,           frac_age,     height,         veget,             &
527                              frac_nobio,    njsc,         veget_max,      fraclut,           &
528                              nwdfraclut,    tot_bare_soil,totfrac_nobio,  qsintmax,          &
529                              co2_flux,      co2_to_bm,    fco2_lu,        temp_growth,       &
530                                                          laydep, rootmass, litter_above, litter_below, carbon, DOC,      &
531                                                          lignin_struc_above,lignin_struc_below, depth_deepsoil)
532    !WRITE(numout,*) 'Sechiba_ini_ZHC2'
533        !WRITE(numout,*) 'litter_above: ', SUM(litter_above(:,:,:,icarbon))
534        !WRITE(numout,*) 'litter_below: ', SUM(litter_below(:,:,:,:,icarbon))
535        !WRITE(numout,*) 'lignin_be',SUM(lignin_struc_below)
536        !WRITE(numout,*) 'DOC_SOC',SUM(DOC),SUM(carbon)
537
538    netco2flux(:) = zero
539    DO jv = 2,nvm
540       netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*veget_max(:,jv)
541    ENDDO
542   
543    !! 1.4 Initialize diffusion coefficients
544    CALL diffuco_initialize (kjit,    kjpindex, index,                  &
545                             rest_id, lalo,     neighbours, resolution, &
546                             rstruct, tq_cdrag)
547   
548    !! 1.5 Initialize remaining variables of energy budget
549        !WRITE(numout,*) 'SUMsechini1_flow',MINVAL(temp_sol(:)),MAXVAL(temp_sol(:)),MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
550    CALL enerbil_initialize (kjit,     kjpindex,     index,    rest_id,  &
551                             qair,                                       &
552                             temp_sol, temp_sol_new, tsol_rad,           &
553                             evapot,   evapot_corr,  qsurf,    fluxsens, &
554                             fluxlat,  vevapp )
555        !WRITE(numout,*) 'SUMsechini2_flow',MINVAL(temp_sol(:)),MAXVAL(temp_sol(:)),MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
556   
557   
558    !! 1.7 Initialize remaining hydrological variables
559    IF ( .NOT. hydrol_cwrr ) THEN
560       ! 1.7.1 Initialize remaining hydrological variables from Choisnel module (2 soil layers)
561       CALL hydrolc_initialize( kjit,      kjpindex,     index,          rest_id, &
562                                veget,     veget_max,    tot_bare_soil,           &
563                                rsol,      drysoil_frac, snow,                    &
564                                snow_age,  snow_nobio,   snow_nobio_age, humrel,  &
565                                vegstress, qsintveg,     shumdiag,       snowrho, &
566                                snowtemp,  snowgrain,    snowdz,         snowheat)
567     
568       evap_bare_lim(:) = -un
569       k_litt(:) = huit
570       
571       ! No specific calculation for shumdiag_perma. We assume it to shumdiag.
572       shumdiag_perma(:,:)=shumdiag(:,:)
573    ELSE
574       !! 1.7.2 Initialize remaining hydrological variables from CWRR module (11 soil layers)
575       CALL hydrol_initialize ( kjit,           kjpindex,  index,         rest_id,          &
576            njsc,           soiltile,  veget,         veget_max,        &
577            humrel,         vegstress, drysoil_frac,                    &
578            shumdiag_perma,    qsintveg,                        &
579            evap_bare_lim,  snow,      snow_age,      snow_nobio,       &
580            snow_nobio_age, snowrho,   snowtemp,      snowgrain,        &
581            snowdz,         snowheat,  &
582            mc_layh,    mcl_layh,  soilmoist)
583
584    ENDIF
585   
586    !! 1.9 Initialize surface parameters (emissivity, albedo and roughness)
587        !WRITE(numout,*) 'SUMsechibaini0',SUM(z0m(:)),SUM(z0h(:)),SUM(zlev(:))
588    CALL condveg_initialize (kjit, kjpindex, index, rest_id, &
589         lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
590         zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
591         drysoil_frac, height, snowdz,snowrho, tot_bare_soil, &
592         temp_air, pb, u, v, lai, &
593         emis, albedo, z0m, z0h, roughheight, &
594         frac_snow_veg,frac_snow_nobio)
595    !WRITE(numout,*) 'SUMsechibaini1',SUM(z0m(:)),SUM(z0h(:)),SUM(zlev(:))
596    !! 1.10 Initialization of soil thermodynamics
597   
598    IF (hydrol_cwrr) THEN
599       CALL thermosoil_initialize (kjit, kjpindex, rest_id,  &
600            temp_sol_new, snow,       shumdiag_perma,        &
601            soilcap,      soilflx,    stempdiag,             &
602            gtemp,               &
603            mc_layh,  mcl_layh,   soilmoist,       njsc ,     &
604            frac_snow_veg,frac_snow_nobio,totfrac_nobio,     &
605            snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb)
606    ELSE
607       CALL thermosoilc_initialize (kjit, kjpindex, rest_id,  &
608            temp_sol_new, snow,       shumdiag_perma,        &
609            soilcap,      soilflx,    stempdiag,             &
610            gtemp, &
611            frac_snow_veg,frac_snow_nobio,totfrac_nobio,     &
612            snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb)
613
614    END IF
615       
616        !! 1.11 Initialize erosion module
617        IF (erosion_module) THEN
618                CALL erosion_initialize(kjit,kjpindex,index,rest_id, lalo,neighbours,resolution, &
619                &       contfrac,std_totsed,gridarea,effgrid_perc)
620        ENDIF !IF (erosion_module) THEN
621    !! 1.12 Initialize river routing
622    IF ( river_routing .AND. nbp_glo .GT. 1) THEN
623       !! 1.12.1 Initialize river routing
624           !WRITE(numout,*) 'SUMsechini3_flow',rest_id,hist_id,hist2_id,MINVAL(riverflow(:,:)),MAXVAL(riverflow(:,:)),MINVAL(coastalflow(:,:)),MAXVAL(coastalflow(:,:))
625       CALL routing_initialize( kjit,        kjpindex,       index,                 &
626            rest_id,     hist_id,        hist2_id,   lalo,      &
627            neighbours,  resolution,     contfrac,   stempdiag, &
628            returnflow,  reinfiltration, irrigation, sed_deposition, poc_deposition, riverflow, &
629            coastalflow, flood_frac,     flood_res, streamfl_frac, stream_frac ,fastr)
630           !WRITE(numout,*) 'SUMsechini4_flow',rest_id,hist_id,hist2_id,MINVAL(riverflow(:,:)),MAXVAL(riverflow(:,:)),MINVAL(coastalflow(:,:)),MAXVAL(coastalflow(:,:))
631    ELSE
632       !! 1.12.2 No routing, set variables to zero
633       riverflow(:,:) = zero
634       coastalflow(:,:) = zero
635       returnflow(:,:) = zero
636       reinfiltration(:,:) = zero
637       irrigation(:,:) = zero
638           sed_deposition(:,:) = zero
639           poc_deposition(:,:) = zero
640       flood_frac(:) = zero
641       streamfl_frac(:) = zero
642       stream_frac(:) = zero
643       flood_res(:) = zero
644       fastr(:) = zero
645    ENDIF
646       
647    !! 1.13 Write internal variables to output fields
648    z0m_out(:) = z0m(:)
649    z0h_out(:) = z0h(:)
650    emis_out(:) = emis(:) 
651    qsurf_out(:) = qsurf(:)
652
653    !! 2. Output variables only once
654    zmaxh_glo(:) = zmaxh
655    CALL xios_orchidee_send_field("zmaxh",zmaxh_glo)
656
657    IF (printlev_loc>=3) WRITE(numout,*) 'sechiba_initialize done'
658
659  END SUBROUTINE sechiba_initialize
660
661!! ==============================================================================================================================\n
662!! SUBROUTINE   : sechiba_main
663!!
664!>\BRIEF        Main routine for the sechiba module performing three functions:
665!! calculating temporal evolution of all variables and preparation of output and
666!! restart files (during the last call only)
667!!
668!!\n DESCRIPTION : Main routine for the sechiba module.
669!! One time step evolution consists of:
670!! - call sechiba_var_init to do some initialization,
671!! - call slowproc_main to do some daily calculations
672!! - call diffuco_main for diffusion coefficient calculation,
673!! - call enerbil_main for energy budget calculation,
674!! - call hydrolc_main (or hydrol_main) for hydrologic processes calculation,
675!! - call enerbil_fusion : last part with fusion,
676!! - call condveg_main for surface conditions such as roughness, albedo, and emmisivity,
677!! - call thermosoil_main(for cwrr) or thermosoilc_main(for choisnel) for soil thermodynamic calculation,
678!! - call sechiba_end to swap previous to new fields.
679!!
680!! RECENT CHANGE(S): None
681!!
682!! MAIN OUTPUT VARIABLE(S): Hydrological variables (:: coastalflow and :: riverflow),
683!! components of the energy budget (:: tsol_rad, :: vevapp, :: fluxsens,
684!! :: temp_sol_new and :: fluxlat), surface characteristics (:: z0_out, :: emis_out,
685!! :: tq_cdrag and :: albedo) and land use related CO2 fluxes (:: netco2flux and
686!! :: fco2_lu)           
687!!
688!! REFERENCE(S) :
689!!
690!! FLOWCHART    :
691!! \latexonly
692!! \includegraphics[scale = 0.5]{sechibamainflow.png}
693!! \endlatexonly
694!! \n
695!_ ================================================================================================================================
696
697  SUBROUTINE sechiba_main (kjit, kjpij, kjpindex, index, &
698       & ldrestart_read, ldrestart_write, &
699       & lalo, contfrac, neighbours, resolution,&
700       & zlev, u, v, qair, q2m, t2m, temp_air, epot_air, ccanopy, &
701       & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
702       & precip_rain, precip_snow, lwdown, swnet, swdown, coszang, pb, &
703       & vevapp, fluxsens, fluxlat, coastalflow, riverflow, netco2flux, fco2_lu, &
704       & tsol_rad, temp_sol_new, qsurf_out, albedo, emis_out, z0m_out, z0h_out,&
705       & veget_out, lai_out, height_out, &
706       & rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
707
708!! 0.1 Input variables
709   
710    INTEGER(i_std), INTENT(in)                               :: kjit              !! Time step number (unitless)
711    INTEGER(i_std), INTENT(in)                               :: kjpij             !! Total size of the un-compressed grid
712                                                                                  !! (unitless)
713    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size - terrestrial pixels only
714                                                                                  !! (unitless)
715    INTEGER(i_std),INTENT (in)                               :: rest_id           !! _Restart_ file identifier (unitless)
716    INTEGER(i_std),INTENT (in)                               :: hist_id           !! _History_ file identifier (unitless)
717    INTEGER(i_std),INTENT (in)                               :: hist2_id          !! _History_ file 2 identifier (unitless)
718    INTEGER(i_std),INTENT (in)                               :: rest_id_stom      !! STOMATE's _Restart_ file identifier
719                                                                                  !! (unitless)
720    INTEGER(i_std),INTENT (in)                               :: hist_id_stom      !! STOMATE's _History_ file identifier
721                                                                                  !! (unitless)
722    INTEGER(i_std),INTENT(in)                                :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file
723                                                                                  !! identifier (unitless)
724    LOGICAL, INTENT(in)                                      :: ldrestart_read    !! Logical for _restart_ file to read
725                                                                                  !! (true/false)
726    LOGICAL, INTENT(in)                                      :: ldrestart_write   !! Logical for _restart_ file to write
727                                                                                  !! (true/false)
728    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)          :: lalo              !! Geographic coordinates (latitude,longitude)
729                                                                                  !! for grid cells (degrees)
730    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: contfrac          !! Fraction of continent in the grid
731                                                                                  !! (unitless, 0-1)
732    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index             !! Indices of the pixels on the map.
733                                                                                  !! Sechiba uses a reduced grid excluding oceans
734                                                                                  !! ::index contains the indices of the
735                                                                                  !! terrestrial pixels only! (unitless)
736    INTEGER(i_std), DIMENSION(kjpindex,NbNeighb), INTENT(in) :: neighbours        !! Neighboring grid points if land!(unitless)
737    REAL(r_std), DIMENSION (kjpindex,2), INTENT(in)          :: resolution        !! Size in x and y of the grid (m)
738   
739    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u                 !! Lowest level wind speed in direction u
740                                                                                  !! @tex $(m.s^{-1})$ @endtex
741    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: v                 !! Lowest level wind speed in direction v
742                                                                                  !! @tex $(m.s^{-1})$ @endtex
743    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: zlev              !! Height of first layer (m)
744    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair              !! Lowest level specific humidity
745                                                                                  !! @tex $(kg kg^{-1})$ @endtex
746    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q2m               !! 2m specific humidity
747                                                                                  !! @tex $(kg kg^{-1})$ @endtex
748    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: t2m               !! 2m air temperature (K)
749    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: precip_rain       !! Rain precipitation
750                                                                                  !! @tex $(kg m^{-2})$ @endtex
751    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: precip_snow       !! Snow precipitation
752                                                                                  !! @tex $(kg m^{-2})$ @endtex
753    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: lwdown            !! Down-welling long-wave flux
754                                                                                  !! @tex $(W m^{-2})$ @endtex
755    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: coszang           !! Cosine of the solar zenith angle (unitless)
756    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swnet             !! Net surface short-wave flux
757                                                                                  !! @tex $(W m^{-2})$ @endtex
758    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swdown            !! Down-welling surface short-wave flux
759                                                                                  !! @tex $(W m^{-2})$ @endtex
760    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_air          !! Air temperature (K)
761    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: epot_air          !! Air potential energy (??J)
762    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ccanopy           !! CO2 concentration in the canopy (ppm)
763    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petAcoef          !! Coefficients A for T from the Planetary
764                                                                                  !! Boundary Layer
765    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqAcoef          !! Coefficients A for q from the Planetary
766                                                                                  !! Boundary Layer
767    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petBcoef          !! Coefficients B for T from the Planetary
768                                                                                  !! Boundary Layer
769    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqBcoef          !! Coefficients B for q from the Planetary
770                                                                                  !! Boundary Layer
771    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: pb                !! Surface pressure (hPa)
772
773
774!! 0.2 Output variables
775
776    REAL(r_std),DIMENSION (kjpindex,nflow), INTENT (out)     :: coastalflow       !! Outflow on coastal points by small basins.
777                                                                                  !! This is the water which flows in a disperse
778                                                                                  !! way into the ocean
779                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
780    REAL(r_std),DIMENSION (kjpindex,nflow), INTENT (out)     :: riverflow         !! Outflow of the major rivers.
781                                                                                  !! The flux will be located on the continental
782                                                                                  !! grid but this should be a coastal point 
783                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
784    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: tsol_rad          !! Radiative surface temperature
785                                                                                  !! @tex $(W m^{-2})$ @endtex
786    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapp            !! Total of evaporation
787                                                                                  !! @tex $(kg m^{-2} days^{-1})$ @endtex
788   
789    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qsurf_out         !! Surface specific humidity
790                                                                                  !! @tex $(kg kg^{-1})$ @endtex
791    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0m_out           !! Surface roughness momentum (output diagnostic, m)
792    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0h_out           !! Surface roughness heat (output diagnostic, m)
793    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out)         :: albedo            !! Surface albedo for visible and near-infrared
794                                                                                  !! (unitless, 0-1)
795    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxsens          !! Sensible heat flux
796                                                                                  !! @tex $(W m^{-2})$ @endtex
797    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxlat           !! Latent heat flux
798                                                                                  !! @tex $(W m^{-2})$ @endtex
799    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: emis_out          !! Emissivity (unitless)
800    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: netco2flux        !! Sum CO2 flux over PFTs
801                                                                                  !! ??(gC m^{-2} s^{-1})??
802    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fco2_lu           !! Land Cover Change CO2 flux
803                                                                                  !! ??(gC m^{-2} s^{-1})??
804    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: veget_out         !! Fraction of vegetation type (unitless, 0-1)
805    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: lai_out           !! Leaf area index (m^2 m^{-2})
806    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: height_out        !! Vegetation Height (m)
807
808!! 0.3 Modified
809
810    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: tq_cdrag          !! Surface drag coefficient  (-)
811    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: temp_sol_new      !! New ground temperature (K)
812
813!! 0.4 local variables
814
815    INTEGER(i_std)                                           :: ji, ig, m,jv, iflow       !! Index (unitless)
816    REAL(r_std), DIMENSION(kjpindex)                         :: histvar           !! Computations for history files (unitless)
817    REAL(r_std), DIMENSION(kjpindex,nlut)                    :: histvar2          !! Computations for history files (unitless)
818    CHARACTER(LEN=80)                                        :: var_name          !! To store variables names for I/O (unitless)
819    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treefrac      !! Total fraction occupied by trees (0-1, uniless)
820    REAL(r_std), DIMENSION(kjpindex)                         :: sum_grassfracC3   !! Total fraction occupied by C3 grasses (0-1, unitless)
821    REAL(r_std), DIMENSION(kjpindex)                         :: sum_grassfracC4   !! Total fraction occupied by C4 grasses (0-1, unitless)
822    REAL(r_std), DIMENSION(kjpindex)                         :: sum_cropfracC3    !! Total fraction occupied by C3 crops (0-1, unitess)
823    REAL(r_std), DIMENSION(kjpindex)                         :: sum_cropfracC4    !! Total fraction occupied by C4 crops (0-1, unitess)
824    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treeFracNdlEvg!! Total fraction occupied by treeFracNdlEvg (0-1, unitess)
825    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treeFracBdlEvg!! Total fraction occupied by treeFracBdlEvg (0-1, unitess)
826    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treeFracNdlDcd!! Total fraction occupied by treeFracNdlDcd (0-1, unitess)
827    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treeFracBdlDcd!! Total fraction occupied by treeFracBdlDcd (0-1, unitess)
828    REAL(r_std), DIMENSION(kjpindex)                         :: grndflux          !! Net energy into soil (W/m2)
829    REAL(r_std), DIMENSION(kjpindex,nsnow)                   :: snowliq           !! Liquid water content (m)
830    REAL(r_std), DIMENSION(kjpindex)                         :: snow_age_diag     !! Only for diag, contains xios_default_val
831    REAL(r_std), DIMENSION(kjpindex,nnobio)                  :: snow_nobio_age_diag !! Only for diag, contains xios_default_val
832    REAL(r_std), DIMENSION(kjpindex)                         :: snowage_glob      !! Snow age on total area including snow on vegetated and bare soil and nobio area @tex ($d$) @endtex
833    REAL(r_std), DIMENSION(kjpindex,nlut)                    :: gpplut            !! GPP on landuse tile, only for diagnostics
834
835
836!_ ================================================================================================================================
837
838    IF (printlev_loc>=3) WRITE(numout,*) 'Start sechiba_main kjpindex =',kjpindex
839
840    !! 1. Initialize variables at each time step
841        !WRITE(numout,*) 'SUMsech0_rau',SUM(rau(:))
842        !WRITE(numout,*) 'SUMsech0_tq_cdrag',SUM(tq_cdrag(:))
843        !WRITE(numout,*) 'SUMsech0_temp_air',SUM(temp_air(:)),temp_air(38)
844        !WRITE(numout,*) 'SUMsech0_contfrac',SUM(contfrac(:)),contfrac(38)
845        !WRITE(numout,*) 'SUMsech0_pb',SUM(pb(:)),pb(38)
846    CALL sechiba_var_init (kjpindex, rau, pb, temp_air) 
847       
848        !WRITE(numout,*) 'SUMsech1_tem_soil',MINVAL(temp_sol(:)),MAXVAL(temp_sol(:)),MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
849        !WRITE(numout,*) 'uv',MINVAL(u(:))      ,MAXVAL(v(:)) ,MINVAL(u(:)),MAXVAL(v(:))
850        !WRITE(numout,*) 'sechiba_check0',rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC
851        !WRITE(numout,*) 'SUMsech1_tq_cdrag',SUM(tq_cdrag(:))
852        !WRITE(numout,*) 'SUMsech1_vu',SUM(v(:)),SUM(u(:))
853    !! 2. Compute diffusion coefficients
854    CALL diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, &
855         & zlev, z0m, z0h, roughheight, temp_sol, temp_air, temp_growth, rau, tq_cdrag, qsurf, qair, q2m, t2m, pb ,  &
856         & rsol, evap_bare_lim, evapot, evapot_corr, snow, flood_frac, flood_res, frac_nobio, snow_nobio, totfrac_nobio, &
857         & swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
858         & vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, gsmean, rveget, rstruct, cimean, gpp, co2_to_bm, &
859         & lalo, neighbours, resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
860         & hist_id, hist2_id)
861    !WRITE(numout,*) 'SUMsech2_tem_soil',MINVAL(temp_sol(:)),MAXVAL(temp_sol(:)),MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
862        !WRITE(numout,*) 'uv',MINVAL(u(:))      ,MAXVAL(v(:)) ,MINVAL(u(:)),MAXVAL(v(:))
863        !WRITE(numout,*) 'sechiba_check1',rest_id, hist_id, hist2_id
864    !! 3. Compute energy balance
865    CALL enerbil_main (kjit, kjpindex, &
866         & index, indexveg, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, &
867         & qair, peqAcoef, peqBcoef, pb, rau, vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
868         & emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
869         & vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo, temp_sol, tsol_rad, &
870         & temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id, &
871         & precip_rain,  pgflux, snowdz, temp_sol_add)
872        !WRITE(numout,*) 'SUMsech3_tem_soil',MINVAL(temp_sol(:))        ,MAXVAL(temp_sol(:)) ,MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
873        !WRITE(numout,*) 'uv',MINVAL(u(:))      ,MAXVAL(v(:)) ,MINVAL(u(:)),MAXVAL(v(:))
874    !! 4. Compute hydrology
875    IF ( .NOT. hydrol_cwrr ) THEN
876       ! 4.1 Water balance from Choisnel module (2 soil layers)
877       CALL hydrolc_main (kjit, kjpindex, index, indexveg, &
878            & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, contfrac, vevapwet, veget, veget_max,&
879            & qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age, tot_melt, transpir, &
880            & precip_rain, precip_snow, returnflow(:,ih2o), reinfiltration(:,ih2o), irrigation(:,ih2o), humrel, vegstress, rsol, drysoil_frac, &
881            & evapot, evapot_corr, flood_frac, flood_res, shumdiag, litterhumdiag, soilcap, rest_id, hist_id, hist2_id, &
882            & temp_air, pb, u, v, swnet, pgflux, &
883            & snowrho,snowtemp,snowgrain,snowdz,snowheat,snowliq, &
884            & grndflux,gtemp,tot_bare_soil, &
885            & lambda_snow,cgrnd_snow,dgrnd_snow,temp_sol_add)
886
887       evap_bare_lim(:) = -un
888       k_litt(:) = huit
889       
890       ! No specific calculation for shumdiag_perma. We assume it to shumdiag.
891       shumdiag_perma(:,:)=shumdiag(:,:)
892    ELSE
893           !WRITE(numout,*) 'Sechiba_floodout0'
894           !WRITE(numout,*) 'floodwater_Min',MINVAL(floodout),MINVAL(runoff),MINVAL(drainage),MINVAL(precip_rain)
895           !WRITE(numout,*) 'floodwater_Max',MAXVAL(floodout),MAXVAL(runoff),MAXVAL(drainage),MAXVAL(precip_rain)
896           !WRITE(numout,*) 'return',MINVAL(returnflow),MAXVAL(returnflow),MINVAL(reinfiltration),MAXVAL(reinfiltration)
897           !DO ig = 1,kjpindex
898                !  IF (floodout(ig).GT.1.0E+10 .OR.floodout(ig).LT.0.0 .OR. ISNAN(floodout(ig))) THEN
899                !     WRITE(numout,*) 'Sechiba_littstr0',ig,floodout(ig)
900                     !WRITE(numout,*) 'Vegetmax',veget_max(ig,:)
901                     !WRITE(numout,*) 'Veget',veget(ig,:)
902                         !WRITE(numout,*) 'totfrac_nobio',totfrac_nobio(ig),'frac_nobio',frac_nobio(ig,:)
903                         !WRITE(numout,*) 'flood_frac',flood_frac(ig)
904                !  ENDIF
905           !ENDDO
906       !! 4.1 Water balance from CWRR module (11 soil layers)
907       CALL hydrol_main (kjit, kjpindex, &
908            & index, indexveg, indexsoil, indexlayer, indexnslm, &
909            & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, &
910            & qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age,  &
911            & tot_melt, transpir, precip_rain, precip_snow, returnflow(:,ih2o), reinfiltration(:,ih2o), irrigation(:,ih2o), &
912            & humrel, vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, flood_frac, flood_res, &
913            & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope,&
914            & rest_id, hist_id, hist2_id,&
915            & contfrac, stempdiag, &
916            & temp_air, pb, u, v, tq_cdrag, swnet, pgflux, &
917            & snowrho,snowtemp,snowgrain,snowdz,snowheat,snowliq, &
918            & grndflux,gtemp,tot_bare_soil, &
919            & lambda_snow,cgrnd_snow,dgrnd_snow,frac_snow_veg,temp_sol_add, &
920            & mc_layh, mcl_layh, soilmoist, runoff_per_soil, drainage_per_soil, soil_mc, &
921            & litter_mc, wat_flux0, wat_flux, precip2canopy, precip2ground, canopy2ground)
922       rsol(:) = -un
923           !WRITE(numout,*) 'SUMsech4_tem_soil',MINVAL(temp_sol(:))     ,MAXVAL(temp_sol(:)) ,MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
924           !WRITE(numout,*) 'Sechiba_floodout1'
925           !WRITE(numout,*) 'floodwater_Min',MINVAL(floodout),MINVAL(runoff),MINVAL(drainage),MINVAL(precip_rain)
926           !WRITE(numout,*) 'floodwater_Max',MAXVAL(floodout),MAXVAL(runoff),MAXVAL(drainage),MAXVAL(precip_rain)
927           !WRITE(numout,*) 'return',MINVAL(returnflow),MAXVAL(returnflow),MINVAL(reinfiltration),MAXVAL(reinfiltration)
928           !DO ig = 1,kjpindex
929                !  IF (floodout(ig).GT.1.0E+10 .OR.floodout(ig).LT.0.0 .OR. ISNAN(floodout(ig))) THEN
930                !     WRITE(numout,*) 'Sechiba_littstr1',ig,floodout(ig)
931                     !WRITE(numout,*) 'Vegetmax',veget_max(ig,:)
932                     !WRITE(numout,*) 'Veget',veget(ig,:)
933                         !WRITE(numout,*) 'totfrac_nobio',totfrac_nobio(ig),'frac_nobio',frac_nobio(ig,:)
934                         !WRITE(numout,*) 'flood_frac',flood_frac(ig)
935                !  ENDIF
936           !ENDDO
937
938    ENDIF
939   
940    !! 5. Compute remaining components of the energy balance
941    IF ( .NOT. ok_explicitsnow ) THEN
942       CALL enerbil_fusion (kjpindex, tot_melt, soilcap, &
943            temp_sol_new, fusion)
944    END IF
945    !WRITE(numout,*) 'SUMsech5_tem_soil',MINVAL(temp_sol(:))    ,MAXVAL(temp_sol(:)) ,MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
946    !! 6. Compute surface variables (emissivity, albedo and roughness)
947    CALL condveg_main (kjit, kjpindex, index, rest_id, hist_id, hist2_id, &
948         lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
949         zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
950         drysoil_frac, height, snowdz, snowrho, tot_bare_soil, &
951         temp_air, pb, u, v, lai, &
952         emis, albedo, z0m, z0h, roughheight, &
953         frac_snow_veg, frac_snow_nobio)
954   
955    !! 7. Compute soil thermodynamics
956    IF (hydrol_cwrr) THEN
957       CALL thermosoil_main (kjit, kjpindex, &
958            index, indexgrnd, &
959            temp_sol_new, snow, soilcap, soilflx, &
960            shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
961            snowdz,snowrho,snowtemp,gtemp,pb,&
962            mc_layh, mcl_layh, soilmoist, njsc,frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
963            lambda_snow, cgrnd_snow, dgrnd_snow)
964                !WRITE(numout,*) 'SUMsech6_tem_soil',MINVAL(temp_sol(:))        ,MAXVAL(temp_sol(:)) ,MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
965    ELSE
966       CALL thermosoilc_main (kjit, kjpindex, &
967            index, indexgrnd, indexnslm, &
968            temp_sol_new, snow, soilcap, soilflx, &
969            shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
970            snowdz,snowrho,snowtemp,gtemp,pb,&
971            frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
972            lambda_snow, cgrnd_snow, dgrnd_snow)
973                !WRITE(numout,*) 'SUMsech7_tem_soil',MINVAL(temp_sol(:))        ,MAXVAL(temp_sol(:)) ,MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
974    END IF
975       
976    !! 8. Compute slow processes (i.e. 'daily' and annual time step)
977    ! ::ok_co2 and ::ok_stomate are flags that determine whether the
978    ! forcing files are written.
979        !WRITE(numout,*) 'Sechiba_ZHC1'
980        !WRITE(numout,*) 'litter_above: ', MAXVAL(litter_above(:,:,:,icarbon))
981        !WRITE(numout,*) 'litter_below: ', MAXVAL(litter_below(:,:,:,:,icarbon))
982        !WRITE(numout,*) 'lignin_be',MAXVAL(lignin_struc_below)
983        !WRITE(numout,*) 'deposition',MAXVAL(sed_deposition),MAXVAL(poc_deposition),MAXVAL(sed_deposition_d),MAXVAL(poc_deposition_d)
984        !WRITE(numout,*) 'DOC_SOC',MAXVAL(carbon),MINVAL(carbon),MAXVAL(DOC),MINVAL(DOC)
985        !WRITE(numout,*) 'floodwater_Min',MINVAL(floodout),MINVAL(runoff),MINVAL(drainage),MINVAL(precip_rain)
986        !WRITE(numout,*) 'floodwater_Max',MAXVAL(floodout),MAXVAL(runoff),MAXVAL(drainage),MAXVAL(precip_rain)
987        !WRITE(numout,*) 'litter_above: ', MAXVAL(litter_above(:,:,:,icarbon))
988        !WRITE(numout,*) 'litter_below: ', MAXVAL(litter_below(:,:,:,:,icarbon))
989        !WRITE(numout,*) 'lignin_be',MAXVAL(lignin_struc_below)
990        !WRITE(numout,*) 'deposition',MAXVAL(sed_deposition),MAXVAL(poc_deposition),MAXVAL(sed_deposition_d),MAXVAL(poc_deposition_d)
991        !WRITE(numout,*) 'DOC_to_topsoil',MAXVAL(DOC_to_topsoil(:,idocl:ipocp)),MAXVAL(DOC_to_subsoil(:,idocl:ipocp)), &
992        !                          MAXVAL(sed_deposition),MAXVAL(poc_deposition), MAXVAL(reinfiltration),MAXVAL(irrigation) 
993        !DO ig = 1,kjpindex
994        !   DO m=1,nvm
995        !       IF (litter_above(ig,2,m,1).GT.1.0E+10 .OR.litter_above(ig,2,m,1).LT.0.0 .OR. ISNAN(litter_above(ig,2,m,1))) THEN
996        !          WRITE(numout,*) 'Sechiba_littstr1',ig,m,litter_above(ig,2,m,1)
997        !       ENDIF
998        !       IF (litter_above(ig,1,m,1).GT.1.0E+10 .OR.litter_above(ig,1,m,1).LT.0.0 .OR. ISNAN(litter_above(ig,1,m,1))) THEN
999        !          WRITE(numout,*) 'Sechiba_littmet1',ig,m,litter_above(ig,1,m,1)
1000        !       ENDIF
1001    !  ENDDO   
1002        !ENDDO
1003        !WRITE(numout,*) 'SUMsech8_tem_soil',MINVAL(temp_sol(:))        ,MAXVAL(temp_sol(:)) ,MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
1004    CALL slowproc_main (kjit, kjpij, kjpindex, &
1005         & index, indexveg, lalo, neighbours, resolution, contfrac, soiltile, fraclut, nwdFraclut, &
1006         & temp_air, temp_sol, stempdiag, textfrac, &
1007         & vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
1008         & deadleaf_cover, &
1009         & assim_param, &
1010         & lai, frac_age, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
1011         & rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
1012         & co2_flux, fco2_lu, co2_to_bm, temp_growth, tot_bare_soil, &
1013         & soil_mc, litter_mc, floodout, runoff, drainage, wat_flux0, wat_flux, &
1014         & drainage_per_soil, runoff_per_soil, DOC_EXP_agg, DOC_to_topsoil, DOC_to_subsoil, &
1015         & (flood_frac+streamfl_frac), stream_frac, precip2canopy, precip2ground, canopy2ground, fastr, &
1016                 & bulkdens,laydep,rootmass,litter_above,litter_below,carbon, DOC,lignin_struc_above, &
1017                 & lignin_struc_below, depth_deepsoil)
1018       
1019           !DO ig = 1,kjpindex
1020                !  IF (floodout(ig).GT.1.0E+10 .OR.floodout(ig).LT.0.0 .OR. ISNAN(floodout(ig))) THEN
1021                !     WRITE(numout,*) 'Sechiba_littstr3',ig,floodout(ig)
1022                     !WRITE(numout,*) 'Vegetmax',veget_max(ig,:)
1023                     !WRITE(numout,*) 'Veget',veget(ig,:)
1024                         !WRITE(numout,*) 'totfrac_nobio',totfrac_nobio(ig),'frac_nobio',frac_nobio(ig,:)
1025                         !WRITE(numout,*) 'flood_frac',flood_frac(ig)
1026                !  ENDIF
1027           !ENDDO
1028        !WRITE(numout,*) 'Sechiba_ZHC2'
1029        !WRITE(numout,*) 'litter_above: ', MAXVAL(litter_above(:,:,:,icarbon))
1030        !WRITE(numout,*) 'litter_below: ', MAXVAL(litter_below(:,:,:,:,icarbon))
1031        !WRITE(numout,*) 'lignin_be',MAXVAL(lignin_struc_below)
1032        !WRITE(numout,*) 'deposition',MAXVAL(sed_deposition),MAXVAL(poc_deposition),MAXVAL(sed_deposition_d),MAXVAL(poc_deposition_d)
1033        !WRITE(numout,*) 'DOC_SOC',MAXVAL(carbon),MINVAL(carbon),MAXVAL(DOC),MINVAL(DOC)
1034        !WRITE(numout,*) 'floodwater_Min',MINVAL(floodout),MINVAL(runoff),MINVAL(drainage),MINVAL(precip_rain)
1035        !WRITE(numout,*) 'floodwater_Max',MAXVAL(floodout),MAXVAL(runoff),MAXVAL(drainage),MAXVAL(precip_rain)
1036        !DO ig = 1,kjpindex
1037        !   DO m=1,nvm
1038        !       IF (litter_above(ig,2,m,1).GT.1.0E+10 .OR.litter_above(ig,2,m,1).LT.0.0 .OR. ISNAN(litter_above(ig,2,m,1))) THEN
1039        !          WRITE(numout,*) 'Sechiba_littstr2',ig,m,litter_above(ig,2,m,1)
1040        !       ENDIF
1041        !       IF (litter_above(ig,1,m,1).GT.1.0E+10 .OR.litter_above(ig,1,m,1).LT.0.0 .OR. ISNAN(litter_above(ig,1,m,1))) THEN
1042        !          WRITE(numout,*) 'Sechiba_littmet2',ig,m,litter_above(ig,1,m,1)
1043        !       ENDIF
1044    !  ENDDO   
1045        !ENDDO
1046        !! Computing soil erosion and updating the vertical profile of soil carbon and litter   ZHC
1047        IF (erosion_module) THEN
1048        !       WRITE(numout,*) 'Erosion_sechiba_main ',carbon(:,1,:,1)
1049            !WRITE(numout,*) 'SUMsech9_tem_soil',MINVAL(temp_sol(:))    ,MAXVAL(temp_sol(:)) ,MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
1050                CALL erosion_main(kjit,kjpindex,rest_id,dt_sechiba,index,ldrestart_read, ldrestart_write, &
1051                        & lalo,neighbours,resolution,temp_sol,bulkdens,textfrac,laydep,contfrac,runoff,drainage,veget, &
1052                        & veget_max,hist_id,std_totsed,gridarea,effgrid_perc,rootmass,litter_above,litter_below,  &
1053                        & carbon, DOC, lignin_struc_above,lignin_struc_below,depth_deepsoil,seddep_rate, &
1054                        & sed_deposition_d,poc_deposition_d,totrunoff_d,peakrunoff_d,totrundra_d,peakrundra_d, &
1055                        & cf_veget,kf_litter,kf_root,erocap,erodep,eroc,tot_erocap,ave_erodep,tot_eroc, &
1056                        & POC_EXP_agg,DOC_ERO_agg,SED_EXP_agg)
1057                !WRITE(numout,*) 'SUMsech10_tem_soil',MINVAL(temp_sol(:))       ,MAXVAL(temp_sol(:)) ,MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
1058!               WRITE(numout,*) 'erosion_main End',erosion_module, sub_erorate(1000000)
1059                !WRITE(numout,*) 'Sechiba_ZHC3'
1060                !WRITE(numout,*) 'litter_above: ', MAXVAL(litter_above(:,:,:,icarbon))
1061                !WRITE(numout,*) 'litter_below: ', MAXVAL(litter_below(:,:,:,:,icarbon))
1062        !       WRITE(numout,*) 'lignin_be',MAXVAL(lignin_struc_below)
1063                !WRITE(numout,*) 'deposition',MAXVAL(sed_deposition),MAXVAL(poc_deposition),MAXVAL(sed_deposition_d),MAXVAL(poc_deposition_d)
1064                !WRITE(numout,*) 'DOC_SOC',MAXVAL(carbon),MINVAL(carbon),MAXVAL(DOC),MINVAL(DOC)
1065                !WRITE(numout,*) 'floodwater_Min',MINVAL(floodout),MINVAL(runoff),MINVAL(drainage),MINVAL(precip_rain)
1066            !WRITE(numout,*) 'floodwater_Max',MAXVAL(floodout),MAXVAL(runoff),MAXVAL(drainage),MAXVAL(precip_rain)
1067        !       WRITE(numout,*) 'floodwater_Max',MAXVAL(totrunoff_d),MAXVAL(peakrunoff_d)
1068        !       DO ig = 1,kjpindex
1069        !      DO m=1,nvm
1070        !           IF (litter_above(ig,2,m,1).GT.1.0E+10 .OR.litter_above(ig,2,m,1).LT.0.0 .OR. ISNAN(litter_above(ig,2,m,1))) THEN
1071        !              WRITE(numout,*) 'Sechiba_littstr2',ig,m,litter_above(ig,2,m,1)
1072        !           ENDIF
1073        !           IF (litter_above(ig,1,m,1).GT.1.0E+10 .OR.litter_above(ig,1,m,1).LT.0.0 .OR. ISNAN(litter_above(ig,1,m,1))) THEN
1074        !              WRITE(numout,*) 'Sechiba_littmet2',ig,m,litter_above(ig,1,m,1)
1075        !           ENDIF
1076    !      ENDDO       
1077        !    ENDDO
1078        ENDIF
1079        !!
1080        !WRITE(numout,*) 'Sechiba_ZHC3'
1081        !DO ig = 1,kjpindex
1082        !       WRITE(numout,*) 'sed_exp',SED_EXP_agg(kjpindex,:)
1083        !       WRITE(numout,*) 'POC_exp',POC_EXP_agg(kjpindex,:)
1084        !       WRITE(numout,*) 'DOC_ero',DOC_ERO_agg(:,iactive)
1085        !       WRITE(numout,*) 'DOC_exp',DOC_EXP_agg(kjpindex,irunoff,:)       
1086        !ENDDO
1087        !WRITE(numout,*) 'SOC',DOC(kjpindex,10,1,1,1,1),carbon(kjpindex,3,10,1)
1088        !WRITE(numout,*) 'sed_exp',SED_EXP_agg(kjpindex,:)
1089        !WRITE(numout,*) 'POC_exp',POC_EXP_agg(kjpindex,:)
1090        !WRITE(numout,*) 'DOC_exp',DOC_EXP_agg(kjpindex,irunoff,:)                                             
1091        !!
1092    !! 9. Compute river routing
1093    IF ( river_routing .AND. nbp_glo .GT. 1) THEN
1094       !! 8.1 River routing
1095          !WRITE(numout,*) 'SUMsech11_tem_soil' ,MINVAL(temp_sol(:))    ,MAXVAL(temp_sol(:)) ,temp_sol(kjpindex)
1096          !WRITE(numout,*) 'Flow11',MINVAL(coastalflow(:,:)),MAXVAL(coastalflow(:,:)),MINVAL(riverflow(:,:)),MAXVAL(riverflow(:,:))
1097          !WRITE(numout,*) 'DOC_SOC',MAXVAL(carbon),MINVAL(carbon),MAXVAL(DOC),MINVAL(DOC)
1098       CALL routing_main (kjit, kjpindex, index, &
1099            & lalo, neighbours, resolution, contfrac, totfrac_nobio, bulkdens, laydep, veget_max, carbon,  &
1100            & depth_deepsoil, floodout, runoff,drainage, transpot, precip_rain, humrel, k_litt, flood_frac,  &
1101            & flood_res, stream_frac, streamfl_frac,stempdiag, reinf_slope, returnflow, reinfiltration,  &
1102                        & irrigation, sed_deposition, poc_deposition,riverflow, coastalflow, rest_id, hist_id, hist2_id, &
1103                        & DOC_EXP_agg, DOC_ERO_agg,POC_EXP_agg, SED_EXP_agg, temp_sol, fastr)
1104            !WRITE(numout,*) 'SUMsech12_tem_soil' ,MINVAL(temp_sol(:)),MAXVAL(temp_sol(:)) ,temp_sol(kjpindex)
1105                !WRITE(numout,*) 'Flow12',MINVAL(coastalflow(:,:)),MAXVAL(coastalflow(:,:)),MINVAL(riverflow(:,:)),MAXVAL(riverflow(:,:))
1106                !WRITE(numout,*) 'DOC_SOC',MAXVAL(carbon),MINVAL(carbon),MAXVAL(DOC),MINVAL(DOC)
1107    ELSE
1108       !! 8.2 No routing, set variables to zero
1109       riverflow(:,:) = zero
1110       coastalflow(:,:) = zero
1111       returnflow(:,:) = zero
1112       reinfiltration(:,:) = zero
1113       irrigation(:,:) = zero
1114           sed_deposition(:,:)= zero
1115           poc_deposition(:,:) = zero
1116       flood_frac(:) = zero
1117       stream_frac(:) = zero
1118       streamfl_frac(:) = zero
1119       flood_res(:) = zero
1120       fastr(:) = zero     
1121!          CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba)
1122!       CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba)
1123    ENDIF
1124   
1125    ! Calculate DOC inputs to soil column from stomate_soilcarbon.f90, to be exported via slowprocesses
1126    DOC_to_topsoil(:,:) = (reinfiltration(:,:)+irrigation(:,:))*mille*one_day/dt_sechiba ! gC/m**2/day
1127    DOC_to_subsoil(:,:) = returnflow(:,:)*mille*one_day/dt_sechiba ! gC/m**2/day
1128        sed_deposition_d(:,:) = sed_deposition(:,:)*one_day/dt_sechiba ! kg/m**2/day
1129        poc_deposition_d(:,:) = poc_deposition(:,:)*one_day/dt_sechiba ! kgC/m**2/day
1130       
1131    !WRITE(numout,*) 'Sechiba_ZHC4'
1132        !WRITE(numout,*) 'flow',MAXVAL(drainage),MAXVAL(transpot)
1133        !WRITE(numout,*) 'deposition',MAXVAL(sed_deposition),MAXVAL(poc_deposition),MAXVAL(sed_deposition_d),MAXVAL(poc_deposition_d)
1134        !WRITE(numout,*) 'floodwater_Min',MAXVAL(reinfiltration(:,idocl:ipocp)),MAXVAL(irrigation(:,idocl:ipocp)),MAXVAL(returnflow(:,idocl:ipocp))
1135        !WRITE(numout,*) 'floodwater_Max',MAXVAL(floodout),MAXVAL(runoff),MAXVAL(drainage),MAXVAL(precip_rain)
1136        !WRITE(numout,*) 'DOC_SOC',MAXVAL(DOC),MAXVAL(carbon),MAXVAL(DOC_to_topsoil),MAXVAL(DOC_to_subsoil)
1137    !! 9.2 Compute global CO2 flux                                             
1138    netco2flux(:) = zero
1139    DO jv = 2,nvm
1140       netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*veget_max(:,jv)
1141    ENDDO
1142 
1143    !! 10. Update the temperature (temp_sol) with newly computed values
1144        !WRITE(numout,*) 'SUMsech13_tem_soil',MINVAL(temp_sol(:)),MAXVAL(temp_sol(:)),MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
1145    CALL sechiba_end (kjpindex, temp_sol_new, temp_sol)
1146        !WRITE(numout,*) 'SUMsech14_tem_soil',MINVAL(temp_sol(:)),MAXVAL(temp_sol(:)),MINVAL(temp_sol_new(:)),MAXVAL(temp_sol_new(:))
1147
1148   
1149    !! 11. Write internal variables to output fields
1150    z0m_out(:) = z0m(:)
1151    z0h_out(:) = z0h(:)
1152    emis_out(:) = emis(:)
1153    qsurf_out(:) = qsurf(:)
1154    veget_out(:,:)  = veget(:,:)
1155    lai_out(:,:)    = lai(:,:)
1156    height_out(:,:) = height(:,:)
1157 
1158    !! 12. Write global variables to history files
1159    sum_treefrac(:) = zero
1160    sum_grassfracC3(:) = zero
1161    sum_grassfracC4(:) = zero
1162    sum_cropfracC3(:) = zero
1163    sum_cropfracC4(:) = zero
1164    sum_treeFracNdlEvg(:) = zero
1165    sum_treeFracBdlEvg(:) = zero
1166    sum_treeFracNdlDcd(:) = zero
1167    sum_treeFracBdlDcd(:) = zero
1168
1169    DO jv = 2, nvm 
1170       IF (is_tree(jv) .AND. natural(jv)) THEN
1171          sum_treefrac(:) = sum_treefrac(:) + veget_max(:,jv)
1172       ELSE IF ((.NOT. is_tree(jv))  .AND. natural(jv)) THEN
1173          ! Grass
1174          IF (is_c4(jv)) THEN
1175             sum_grassfracC4(:) = sum_grassfracC4(:) + veget_max(:,jv)
1176          ELSE
1177             sum_grassfracC3(:) = sum_grassfracC3(:) + veget_max(:,jv)
1178          END IF
1179       ELSE 
1180          ! Crop and trees not natural
1181          IF (is_c4(jv)) THEN
1182             sum_cropfracC4(:) = sum_cropfracC4(:) + veget_max(:,jv)
1183          ELSE
1184             sum_cropfracC3(:) = sum_cropfracC3(:) + veget_max(:,jv)
1185          END IF
1186       ENDIF
1187
1188       IF (is_tree(jv)) THEN
1189          IF (is_evergreen(jv)) THEN
1190             IF (is_needleleaf(jv)) THEN
1191                ! Fraction for needleleaf evergreen trees (treeFracNdlEvg)
1192                sum_treeFracNdlEvg(:) = sum_treeFracNdlEvg(:) + veget_max(:,jv)
1193             ELSE
1194                ! Fraction for broadleaf evergreen trees (treeFracBdlEvg)
1195                sum_treeFracBdlEvg(:) = sum_treeFracBdlEvg(:) + veget_max(:,jv)
1196             END IF
1197          ELSE IF (is_deciduous(jv)) THEN
1198             IF (is_needleleaf(jv)) THEN
1199                ! Fraction for needleleaf deciduous trees (treeFracNdlDcd)
1200                sum_treeFracNdlDcd(:) = sum_treeFracNdlDcd(:) + veget_max(:,jv)
1201             ELSE 
1202                ! Fraction for broadleafs deciduous trees (treeFracBdlDcd)
1203                sum_treeFracBdlDcd(:) = sum_treeFracBdlDcd(:) + veget_max(:,jv)
1204             END IF
1205          END IF
1206       END IF
1207    ENDDO         
1208
1209    histvar(:)=zero
1210    DO jv = 2, nvm
1211       IF (is_deciduous(jv)) THEN
1212          histvar(:) = histvar(:) + veget_max(:,jv)*100*contfrac
1213       ENDIF
1214    ENDDO
1215    CALL xios_orchidee_send_field("treeFracPrimDec",histvar)
1216
1217    histvar(:)=zero
1218    DO jv = 2, nvm
1219       IF (is_evergreen(jv)) THEN
1220          histvar(:) = histvar(:) + veget_max(:,jv)*100*contfrac
1221       ENDIF
1222    ENDDO
1223    CALL xios_orchidee_send_field("treeFracPrimEver",histvar)
1224
1225    histvar(:)=zero
1226    DO jv = 2, nvm
1227       IF ( .NOT.(is_c4(jv)) ) THEN
1228          histvar(:) = histvar(:) + veget_max(:,jv)*100*contfrac
1229       ENDIF
1230    ENDDO
1231    CALL xios_orchidee_send_field("c3PftFrac",histvar)
1232
1233    histvar(:)=zero
1234    DO jv = 2, nvm
1235       IF ( is_c4(jv) ) THEN
1236          histvar(:) = histvar(:) + veget_max(:,jv)*100*contfrac
1237       ENDIF
1238    ENDDO
1239    CALL xios_orchidee_send_field("c4PftFrac",histvar)
1240
1241    CALL xios_orchidee_send_field("temp_sol_new",temp_sol_new)
1242    CALL xios_orchidee_send_field("fluxsens",fluxsens)
1243    CALL xios_orchidee_send_field("fluxlat",fluxlat)
1244
1245
1246    ! Add XIOS default value where no snow
1247    DO ji=1,kjpindex 
1248       IF (snow(ji) .GT. zero) THEN
1249          snow_age_diag(ji) = snow_age(ji)
1250          snow_nobio_age_diag(ji,:) = snow_nobio_age(ji,:)
1251       
1252          snowage_glob(ji) = snow_age(ji)*frac_snow_veg(ji)*(1-totfrac_nobio(ji)) + &
1253               SUM(snow_nobio_age(ji,:)*frac_snow_nobio(ji,:)*frac_nobio(ji,:))
1254          IF (snowage_glob(ji) .NE. 0) snowage_glob(ji) = snowage_glob(ji) / &
1255               (frac_snow_veg(ji)*(1-totfrac_nobio(ji)) + SUM(frac_snow_nobio(ji,:)*frac_nobio(ji,:)))
1256       ELSE
1257          snow_age_diag(ji) = xios_default_val
1258          snow_nobio_age_diag(ji,:) = xios_default_val
1259          snowage_glob(ji) = xios_default_val
1260       END IF
1261    END DO
1262   
1263    CALL xios_orchidee_send_field("snow",snow)
1264    CALL xios_orchidee_send_field("snowage",snow_age_diag)
1265    CALL xios_orchidee_send_field("flood_frac",flood_frac)
1266    CALL xios_orchidee_send_field("stream_frac",stream_frac)
1267    CALL xios_orchidee_send_field("streamfl_frac",streamfl_frac)                               
1268    CALL xios_orchidee_send_field("snownobio",snow_nobio)
1269    CALL xios_orchidee_send_field("snownobioage",snow_nobio_age_diag)
1270    CALL xios_orchidee_send_field("snowage_glob",snowage_glob)
1271
1272    CALL xios_orchidee_send_field("frac_snow", SUM(frac_snow_nobio,2)*totfrac_nobio+frac_snow_veg*(1-totfrac_nobio))
1273    CALL xios_orchidee_send_field("frac_snow_veg", frac_snow_veg)
1274    CALL xios_orchidee_send_field("frac_snow_nobio", frac_snow_nobio)
1275    CALL xios_orchidee_send_field("reinf_slope",reinf_slope)
1276    CALL xios_orchidee_send_field("njsc",REAL(njsc, r_std))
1277    CALL xios_orchidee_send_field("vegetfrac",veget)
1278    CALL xios_orchidee_send_field("maxvegetfrac",veget_max)
1279    CALL xios_orchidee_send_field("nobiofrac",frac_nobio)
1280    CALL xios_orchidee_send_field("soiltile",soiltile)
1281    CALL xios_orchidee_send_field("rstruct",rstruct)
1282    CALL xios_orchidee_send_field("gpp",gpp/dt_sechiba)
1283    CALL xios_orchidee_send_field("gpp_ipcc2",SUM(gpp,dim=2)/dt_sechiba)
1284
1285    histvar(:)=zero
1286    DO jv = 2, nvm
1287       IF ( .NOT. is_tree(jv) .AND. natural(jv) ) THEN
1288          histvar(:) = histvar(:) + gpp(:,jv)
1289       ENDIF
1290    ENDDO
1291    CALL xios_orchidee_send_field("gppgrass",histvar/dt_sechiba)
1292
1293    histvar(:)=zero
1294    DO jv = 2, nvm
1295       IF ( (.NOT. is_tree(jv)) .AND. (.NOT. natural(jv)) ) THEN
1296          histvar(:) = histvar(:) + gpp(:,jv)
1297       ENDIF
1298    ENDDO
1299    CALL xios_orchidee_send_field("gppcrop",histvar/dt_sechiba)
1300
1301    histvar(:)=zero
1302    DO jv = 2, nvm
1303       IF ( is_tree(jv) ) THEN
1304          histvar(:) = histvar(:) + gpp(:,jv)
1305       ENDIF
1306    ENDDO
1307    CALL xios_orchidee_send_field("gpptree",histvar/dt_sechiba)
1308
1309    CALL xios_orchidee_send_field("nee",co2_flux/dt_sechiba)
1310    CALL xios_orchidee_send_field("drysoil_frac",drysoil_frac)
1311    CALL xios_orchidee_send_field("vevapflo",vevapflo/dt_sechiba)
1312    CALL xios_orchidee_send_field("k_litt",k_litt)
1313    CALL xios_orchidee_send_field("beta",vbeta)
1314    CALL xios_orchidee_send_field("vbeta1",vbeta1)
1315    CALL xios_orchidee_send_field("vbeta2",vbeta2)
1316    CALL xios_orchidee_send_field("vbeta3",vbeta3)
1317    CALL xios_orchidee_send_field("vbeta4",vbeta4)
1318    CALL xios_orchidee_send_field("vbeta5",vbeta5)
1319    CALL xios_orchidee_send_field("gsmean",gsmean)
1320    CALL xios_orchidee_send_field("cimean",cimean)
1321    CALL xios_orchidee_send_field("rveget",rveget)
1322    CALL xios_orchidee_send_field("rsol",rsol)
1323 
1324    histvar(:)=SUM(vevapwet(:,:),dim=2)
1325    CALL xios_orchidee_send_field("evspsblveg",histvar/dt_sechiba)
1326    histvar(:)= vevapnu(:)+vevapsno(:)
1327    CALL xios_orchidee_send_field("evspsblsoi",histvar/dt_sechiba)
1328    histvar(:)=SUM(transpir(:,:),dim=2)
1329    CALL xios_orchidee_send_field("tran",histvar/dt_sechiba)
1330
1331    ! For CMIP6 data request: the following fractions are fractions of the total grid-cell,
1332    ! which explains the multiplication by contfrac
1333    CALL xios_orchidee_send_field("treeFrac",sum_treefrac(:)*100*contfrac(:))
1334    CALL xios_orchidee_send_field("grassFracC3",sum_grassfracC3(:)*100*contfrac(:))
1335    CALL xios_orchidee_send_field("grassFracC4",sum_grassfracC4(:)*100*contfrac(:))
1336    CALL xios_orchidee_send_field("cropFracC3",sum_cropfracC3(:)*100*contfrac(:))
1337    CALL xios_orchidee_send_field("cropFracC4",sum_cropfracC4(:)*100*contfrac(:))
1338    CALL xios_orchidee_send_field("treeFracNdlEvg",sum_treeFracNdlEvg(:)*100*contfrac(:))
1339    CALL xios_orchidee_send_field("treeFracBdlEvg",sum_treeFracBdlEvg(:)*100*contfrac(:))
1340    CALL xios_orchidee_send_field("treeFracNdlDcd",sum_treeFracNdlDcd(:)*100*contfrac(:))
1341    CALL xios_orchidee_send_field("treeFracBdlDcd",sum_treeFracBdlDcd(:)*100*contfrac(:))
1342
1343    histvar(:)=veget_max(:,1)*100*contfrac(:)
1344    CALL xios_orchidee_send_field("baresoilFrac",histvar)
1345    histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)*100*contfrac(:)
1346    CALL xios_orchidee_send_field("residualFrac",histvar)
1347
1348    ! For CMIP6 data request: cnc = canopy cover fraction over land area
1349    histvar(:)=zero
1350    DO jv=2,nvm
1351       histvar(:) = histvar(:) + veget_max(:,jv)*100
1352    END DO
1353    CALL xios_orchidee_send_field("cnc",histvar)
1354   
1355    CALL xios_orchidee_send_field("tsol_rad",tsol_rad-273.15)
1356    CALL xios_orchidee_send_field("qsurf",qsurf)
1357    CALL xios_orchidee_send_field("emis",emis)
1358    CALL xios_orchidee_send_field("z0m",z0m)
1359    CALL xios_orchidee_send_field("z0h",z0h)
1360    CALL xios_orchidee_send_field("roughheight",roughheight)
1361    CALL xios_orchidee_send_field("lai",lai)
1362    histvar(:)=zero   
1363    DO ji = 1, kjpindex
1364       IF (SUM(veget_max(ji,:)) > zero) THEN
1365         DO jv=2,nvm
1366            histvar(ji) = histvar(ji) + veget_max(ji,jv)*lai(ji,jv)/SUM(veget_max(ji,:))
1367         END DO
1368       END IF
1369    END DO
1370    CALL xios_orchidee_send_field("LAImean",histvar)
1371    CALL xios_orchidee_send_field("vevapsno",vevapsno/dt_sechiba)
1372    CALL xios_orchidee_send_field("vevapp",vevapp/dt_sechiba)
1373    CALL xios_orchidee_send_field("vevapnu",vevapnu/dt_sechiba)
1374    CALL xios_orchidee_send_field("transpir",transpir*one_day/dt_sechiba)
1375    CALL xios_orchidee_send_field("inter",vevapwet*one_day/dt_sechiba)
1376    CALL xios_orchidee_send_field("Qf",fusion)
1377    histvar(:)=zero
1378    DO jv=1,nvm
1379      histvar(:) = histvar(:) + vevapwet(:,jv)
1380    ENDDO
1381    CALL xios_orchidee_send_field("ECanop",histvar/dt_sechiba)
1382    histvar(:)=zero
1383    DO jv=1,nvm
1384      histvar(:) = histvar(:) + transpir(:,jv)
1385    ENDDO
1386    CALL xios_orchidee_send_field("TVeg",histvar/dt_sechiba)
1387    CALL xios_orchidee_send_field("ACond",tq_cdrag)
1388
1389
1390    !! Calculate diagnostic variables on Landuse tiles for LUMIP/CMIP6
1391
1392    ! Calculate fraction of landuse tiles related to the whole grid cell
1393    DO jv=1,nlut
1394       histvar2(:,jv) = fraclut(:,jv) * contfrac(:)
1395    END DO
1396    CALL xios_orchidee_send_field("fraclut",histvar2)
1397
1398    CALL xios_orchidee_send_field("nwdFraclut",nwdFraclut(:,:))
1399   
1400    ! Calculate GPP on landuse tiles
1401    ! val_exp is used as missing value where the values are not known i.e. where the tile is not represented
1402    ! or for pasture (id_pst) or urban land (id_urb).
1403    gpplut(:,:)=0
1404    DO jv=1,nvm
1405       IF (natural(jv)) THEN
1406          gpplut(:,id_psl) = gpplut(:,id_psl) + gpp(:,jv)
1407       ELSE
1408          gpplut(:,id_crp) = gpplut(:,id_crp) + gpp(:,jv)
1409       ENDIF
1410    END DO
1411
1412    ! Transform from gC/m2/s into kgC/m2/s
1413    WHERE (fraclut(:,id_psl)>min_sechiba)
1414       gpplut(:,id_psl) = gpplut(:,id_psl)/fraclut(:,id_psl)/1000
1415    ELSEWHERE
1416       gpplut(:,id_psl) = xios_default_val
1417    END WHERE
1418    WHERE (fraclut(:,id_crp)>min_sechiba)
1419       gpplut(:,id_crp) = gpplut(:,id_crp)/fraclut(:,id_crp)/1000
1420    ELSEWHERE
1421       gpplut(:,id_crp) = xios_default_val
1422    END WHERE
1423    gpplut(:,id_pst) = xios_default_val
1424    gpplut(:,id_urb) = xios_default_val
1425
1426    CALL xios_orchidee_send_field("gpplut",gpplut)
1427
1428
1429    IF ( .NOT. almaoutput ) THEN
1430       ! Write history file in IPSL-format
1431       CALL histwrite_p(hist_id, 'beta', kjit, vbeta, kjpindex, index)
1432       CALL histwrite_p(hist_id, 'z0m', kjit, z0m, kjpindex, index)
1433       CALL histwrite_p(hist_id, 'z0h', kjit, z0h, kjpindex, index)
1434       CALL histwrite_p(hist_id, 'roughheight', kjit, roughheight, kjpindex, index)
1435       CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1436       CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1437       CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1438       CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1439       CALL histwrite_p(hist_id, 'subli', kjit, vevapsno, kjpindex, index)
1440       CALL histwrite_p(hist_id, 'evapnu', kjit, vevapnu, kjpindex, index)
1441       CALL histwrite_p(hist_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
1442       CALL histwrite_p(hist_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
1443       CALL histwrite_p(hist_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
1444       CALL histwrite_p(hist_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
1445       CALL histwrite_p(hist_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
1446       CALL histwrite_p(hist_id, 'vbeta4', kjit, vbeta4, kjpindex, index)   
1447       CALL histwrite_p(hist_id, 'vbeta5', kjit, vbeta5, kjpindex, index)   
1448       CALL histwrite_p(hist_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
1449       CALL histwrite_p(hist_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
1450       CALL histwrite_p(hist_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
1451
1452       IF ( .NOT. hydrol_cwrr ) THEN
1453          CALL histwrite_p(hist_id, 'rsol', kjit, rsol, kjpindex, index)
1454       ENDIF
1455       CALL histwrite_p(hist_id, 'snow', kjit, snow, kjpindex, index)
1456       CALL histwrite_p(hist_id, 'snowage', kjit, snow_age, kjpindex, index)
1457       CALL histwrite_p(hist_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
1458       CALL histwrite_p(hist_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
1459
1460       IF (ok_explicitsnow) THEN
1461          CALL histwrite_p(hist_id, 'grndflux', kjit, grndflux, kjpindex,index)
1462          CALL histwrite_p(hist_id, 'snowtemp',kjit,snowtemp,kjpindex*nsnow,indexsnow)
1463          CALL histwrite_p(hist_id, 'snowliq', kjit,snowliq,kjpindex*nsnow,indexsnow)
1464          CALL histwrite_p(hist_id, 'snowdz', kjit,snowdz,kjpindex*nsnow,indexsnow)
1465          CALL histwrite_p(hist_id, 'snowrho', kjit,snowrho,kjpindex*nsnow,indexsnow)
1466          CALL histwrite_p(hist_id, 'snowgrain',kjit,snowgrain,kjpindex*nsnow,indexsnow)
1467          CALL histwrite_p(hist_id, 'snowheat',kjit,snowheat,kjpindex*nsnow,indexsnow)
1468       END IF
1469
1470       CALL histwrite_p(hist_id,'frac_snow_veg',kjit,frac_snow_veg,kjpindex,index)
1471       CALL histwrite_p(hist_id, 'frac_snow_nobio', kjit,frac_snow_nobio,kjpindex*nnobio, indexnobio)
1472       CALL histwrite_p(hist_id, 'pgflux',kjit,pgflux,kjpindex,index)
1473       CALL histwrite_p(hist_id, 'soiltile',  kjit, soiltile, kjpindex*nstm, indexsoil)
1474       !
1475       IF ( hydrol_cwrr ) THEN
1476          CALL histwrite_p(hist_id, 'soilindex',  kjit, REAL(njsc, r_std), kjpindex, index)
1477          CALL histwrite_p(hist_id, 'reinf_slope',  kjit, reinf_slope, kjpindex, index)
1478          CALL histwrite_p(hist_id, 'k_litt', kjit, k_litt, kjpindex, index)
1479       ENDIF
1480       IF ( do_floodplains ) THEN
1481          CALL histwrite_p(hist_id, 'evapflo', kjit, vevapflo, kjpindex, index)
1482          CALL histwrite_p(hist_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
1483       ENDIF
1484       IF ( ok_co2 ) THEN
1485          CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)   
1486          CALL histwrite_p(hist_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
1487          CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)   
1488       ENDIF
1489       IF ( ok_stomate ) THEN
1490          CALL histwrite_p(hist_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)   
1491       ENDIF
1492
1493       histvar(:)=SUM(vevapwet(:,:),dim=2)
1494       CALL histwrite_p(hist_id, 'evspsblveg', kjit, histvar, kjpindex, index)
1495
1496       histvar(:)= vevapnu(:)+vevapsno(:)
1497       CALL histwrite_p(hist_id, 'evspsblsoi', kjit, histvar, kjpindex, index)
1498
1499       histvar(:)=SUM(transpir(:,:),dim=2)
1500       CALL histwrite_p(hist_id, 'tran', kjit, histvar, kjpindex, index)
1501
1502       histvar(:)= sum_treefrac(:)*100*contfrac(:)
1503       CALL histwrite_p(hist_id, 'treeFrac', kjit, histvar, kjpindex, index) 
1504
1505       histvar(:)= (sum_grassfracC3(:)+sum_grassfracC4(:))*100*contfrac(:)
1506       CALL histwrite_p(hist_id, 'grassFrac', kjit, histvar, kjpindex, index) 
1507
1508       histvar(:)= (sum_cropfracC3(:)+sum_cropfracC4(:))*100*contfrac(:)
1509       CALL histwrite_p(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)
1510
1511       histvar(:)=veget_max(:,1)*100*contfrac(:)
1512       CALL histwrite_p(hist_id, 'baresoilFrac', kjit, histvar, kjpindex, index)
1513
1514       histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)*100*contfrac(:)
1515       CALL histwrite_p(hist_id, 'residualFrac', kjit, histvar, kjpindex, index)
1516    ELSE
1517       ! Write history file in ALMA format
1518       CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1519       CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1520       CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1521       CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1522       CALL histwrite_p(hist_id, 'Qf', kjit, fusion, kjpindex, index)
1523       CALL histwrite_p(hist_id, 'ESoil', kjit, vevapnu, kjpindex, index)
1524       CALL histwrite_p(hist_id, 'EWater', kjit, vevapflo, kjpindex, index)
1525       CALL histwrite_p(hist_id, 'SWE', kjit, snow, kjpindex, index)
1526       histvar(:)=zero
1527       DO jv=1,nvm
1528          histvar(:) = histvar(:) + transpir(:,jv)
1529       ENDDO
1530       CALL histwrite_p(hist_id, 'TVeg', kjit, histvar, kjpindex, index)
1531       histvar(:)=zero
1532       DO jv=1,nvm
1533          histvar(:) = histvar(:) + vevapwet(:,jv)
1534       ENDDO
1535       CALL histwrite_p(hist_id, 'ECanop', kjit, histvar, kjpindex, index)
1536       CALL histwrite_p(hist_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
1537       CALL histwrite_p(hist_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
1538       !
1539       CALL histwrite_p(hist_id, 'Z0m', kjit, z0m, kjpindex, index)
1540       CALL histwrite_p(hist_id, 'Z0h', kjit, z0h, kjpindex, index)
1541       CALL histwrite_p(hist_id, 'EffectHeight', kjit, roughheight, kjpindex, index)
1542       !
1543       IF ( do_floodplains ) THEN
1544          CALL histwrite_p(hist_id, 'Qflood', kjit, vevapflo, kjpindex, index)
1545          CALL histwrite_p(hist_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
1546       ENDIF
1547       !
1548       IF ( ok_co2 ) THEN
1549          CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)   
1550          CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)   
1551          CALL histwrite_p(hist_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
1552       ENDIF
1553       IF ( ok_stomate ) THEN
1554             CALL histwrite_p(hist_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)   
1555       ENDIF
1556    ENDIF ! almaoutput
1557   
1558    !! 13. Write additional output file with higher frequency
1559    IF ( hist2_id > 0 ) THEN
1560       IF ( .NOT. almaoutput ) THEN
1561          ! Write history file in IPSL-format
1562          CALL histwrite_p(hist2_id, 'tsol_rad', kjit, tsol_rad, kjpindex, index)
1563          CALL histwrite_p(hist2_id, 'qsurf', kjit, qsurf, kjpindex, index)
1564          CALL histwrite_p(hist2_id, 'albedo', kjit, albedo, kjpindex*2, indexalb)
1565          CALL histwrite_p(hist2_id, 'emis', kjit, emis, kjpindex, index)
1566          CALL histwrite_p(hist2_id, 'beta', kjit, vbeta, kjpindex, index)
1567          CALL histwrite_p(hist2_id, 'z0m', kjit, z0m, kjpindex, index)
1568          CALL histwrite_p(hist2_id, 'z0h', kjit, z0h, kjpindex, index)
1569          CALL histwrite_p(hist2_id, 'roughheight', kjit, roughheight, kjpindex, index)
1570          CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1571          CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1572          CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1573          CALL histwrite_p(hist2_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1574          CALL histwrite_p(hist2_id, 'subli', kjit, vevapsno, kjpindex, index)
1575          IF ( do_floodplains ) THEN
1576             CALL histwrite_p(hist2_id, 'vevapflo', kjit, vevapflo, kjpindex, index)
1577             CALL histwrite_p(hist2_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
1578          ENDIF
1579          CALL histwrite_p(hist2_id, 'vevapnu', kjit, vevapnu, kjpindex, index)
1580          CALL histwrite_p(hist2_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
1581          CALL histwrite_p(hist2_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
1582          CALL histwrite_p(hist2_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
1583          CALL histwrite_p(hist2_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
1584          CALL histwrite_p(hist2_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
1585          CALL histwrite_p(hist2_id, 'vbeta4', kjit, vbeta4, kjpindex, index)   
1586          CALL histwrite_p(hist2_id, 'vbeta5', kjit, vbeta5, kjpindex, index)   
1587          CALL histwrite_p(hist2_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
1588          CALL histwrite_p(hist2_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
1589          CALL histwrite_p(hist2_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
1590          IF ( .NOT. hydrol_cwrr ) THEN
1591             CALL histwrite_p(hist2_id, 'rsol', kjit, rsol, kjpindex, index)
1592          ENDIF
1593          CALL histwrite_p(hist2_id, 'snow', kjit, snow, kjpindex, index)
1594          CALL histwrite_p(hist2_id, 'snowage', kjit, snow_age, kjpindex, index)
1595          CALL histwrite_p(hist2_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
1596          CALL histwrite_p(hist2_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
1597          !
1598          IF (  hydrol_cwrr ) THEN
1599             CALL histwrite_p(hist2_id, 'soilindex',  kjit, REAL(njsc, r_std), kjpindex, index)
1600             CALL histwrite_p(hist2_id, 'reinf_slope',  kjit, reinf_slope, kjpindex, index)
1601          ENDIF
1602          !
1603          IF ( ok_co2 ) THEN
1604             CALL histwrite_p(hist2_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)   
1605             CALL histwrite_p(hist2_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
1606             CALL histwrite_p(hist2_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)   
1607          ENDIF
1608          IF ( ok_stomate ) THEN
1609             CALL histwrite_p(hist2_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)   
1610          ENDIF
1611       ELSE
1612          ! Write history file in ALMA format
1613          CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1614          CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1615          CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1616          CALL histwrite_p(hist2_id, 'Qf', kjit, fusion, kjpindex, index)
1617          CALL histwrite_p(hist2_id, 'ESoil', kjit, vevapnu, kjpindex, index)
1618          IF ( do_floodplains ) THEN
1619             CALL histwrite_p(hist2_id, 'EWater', kjit, vevapflo, kjpindex, index)
1620             CALL histwrite_p(hist2_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
1621          ENDIF
1622          CALL histwrite_p(hist2_id, 'SWE', kjit, snow, kjpindex, index)
1623          histvar(:)=zero
1624          DO jv=1,nvm
1625             histvar(:) = histvar(:) + transpir(:,jv)
1626          ENDDO
1627          CALL histwrite_p(hist2_id, 'TVeg', kjit, histvar, kjpindex, index)
1628          histvar(:)=zero
1629          DO jv=1,nvm
1630             histvar(:) = histvar(:) + vevapwet(:,jv)
1631          ENDDO
1632          CALL histwrite_p(hist2_id, 'ECanop', kjit, histvar, kjpindex, index)
1633          CALL histwrite_p(hist2_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
1634          CALL histwrite_p(hist2_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
1635          IF ( ok_co2 ) THEN
1636             CALL histwrite_p(hist2_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
1637          ENDIF
1638          IF ( ok_stomate ) THEN
1639             CALL histwrite_p(hist2_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)   
1640          ENDIF
1641       ENDIF ! almaoutput
1642    ENDIF ! hist2_id
1643
1644
1645    !! Change the vegetation fractions if a new map was read in slowproc. This is done
1646    !! after lcchange has been done in stomatelpj
1647    IF (done_stomate_lcchange) THEN
1648       CALL slowproc_change_frac(kjpindex, lai, &
1649                                 veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile, fraclut, nwdFraclut)
1650       done_stomate_lcchange=.FALSE.
1651    END IF
1652
1653    !! 14. If it is the last time step, write restart files
1654    IF (ldrestart_write) THEN
1655       CALL sechiba_finalize( &
1656            kjit,     kjpij,  kjpindex, index,   rest_id, &
1657            tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad)
1658    END IF
1659
1660  END SUBROUTINE sechiba_main
1661
1662
1663!!  =============================================================================================================================
1664!! SUBROUTINE:    sechiba_finalize
1665!!
1666!>\BRIEF          Finalize all modules by calling their "_finalize" subroutines.
1667!!
1668!! DESCRIPTION:   Finalize all modules by calling their "_finalize" subroutines. These subroutines will write variables to
1669!!                restart file.
1670!!
1671!! \n
1672!_ ==============================================================================================================================
1673
1674  SUBROUTINE sechiba_finalize( &
1675       kjit,     kjpij,  kjpindex, index,   rest_id, &
1676       tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad)
1677
1678!! 0.1 Input variables   
1679    INTEGER(i_std), INTENT(in)                               :: kjit              !! Time step number (unitless)
1680    INTEGER(i_std), INTENT(in)                               :: kjpij             !! Total size of the un-compressed grid
1681                                                                                  !! (unitless)
1682    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size - terrestrial pixels only
1683                                                                                  !! (unitless)
1684    INTEGER(i_std),INTENT (in)                               :: rest_id           !! _Restart_ file identifier (unitless)
1685    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index             !! Indices of the pixels on the map.
1686                                                                                  !! Sechiba uses a reduced grid excluding oceans
1687                                                                                  !! ::index contains the indices of the
1688                                                                                  !! terrestrial pixels only! (unitless)
1689    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: tsol_rad           !! Radiative surface temperature
1690                                                                                  !! @tex $(W m^{-2})$ @endtex
1691    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: vevapp             !! Total of evaporation
1692                                                                                  !! @tex $(kg m^{-2} days^{-1})$ @endtex
1693    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: fluxsens           !! Sensible heat flux
1694                                                                                  !! @tex $(W m^{-2})$ @endtex
1695    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: fluxlat            !! Latent heat flux
1696                                                                                  !! @tex $(W m^{-2})$ @endtex
1697    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: tq_cdrag           !! Surface drag coefficient (-)
1698
1699!! 0.2 Local variables
1700    INTEGER(i_std)                                          :: ji, jv             !! Index (unitless)
1701    REAL(r_std), DIMENSION(kjpindex)                        :: histvar            !! Computations for history files (unitless)
1702    CHARACTER(LEN=80)                                       :: var_name           !! To store variables names for I/O (unitless)
1703
1704
1705    !! Write restart file for the next simulation from SECHIBA and other modules
1706
1707    IF (printlev_loc>=3) WRITE (numout,*) 'Start sechiba_finalize for writing restart files'
1708
1709    !! 1. Call diffuco_finalize to write restart files
1710    CALL diffuco_finalize (kjit, kjpindex, rest_id, rstruct )
1711   
1712    !! 2. Call energy budget to write restart files
1713    CALL enerbil_finalize (kjit,   kjpindex,    rest_id,            &
1714                           evapot, evapot_corr, temp_sol, tsol_rad, &
1715                           qsurf,  fluxsens,    fluxlat,  vevapp )
1716   
1717    !! 3. Call hydrology to write restart files
1718    IF ( .NOT. hydrol_cwrr ) THEN
1719       !! 3.1 Call water balance from Choisnel module (2 soil layers) to write restart file
1720       CALL hydrolc_finalize( kjit,      kjpindex,   rest_id,        snow,       &
1721                              snow_age,  snow_nobio, snow_nobio_age, humrel,     &
1722                              vegstress, qsintveg,   snowrho,        snowtemp,   &
1723                              snowdz,    snowheat,   snowgrain,       &
1724                              drysoil_frac, rsol, shumdiag)
1725
1726       evap_bare_lim(:) = -un
1727       k_litt(:) = huit
1728       shumdiag_perma(:,:)=shumdiag(:,:)
1729    ELSE
1730       !! 3.2 Call water balance from CWRR module (11 soil layers) to write restart file
1731       CALL hydrol_finalize( kjit,           kjpindex, rest_id,  vegstress,  &
1732                             qsintveg,       humrel,   snow,     snow_age, snow_nobio, &
1733                             snow_nobio_age, snowrho,  snowtemp, snowdz,     &
1734                             snowheat,       snowgrain,    &
1735                             drysoil_frac, evap_bare_lim)
1736    ENDIF
1737   
1738    !! 4. Call condveg to write surface variables to restart files
1739    CALL condveg_finalize (kjit, kjpindex, rest_id, z0m, z0h, roughheight)
1740   
1741    !! 5. Call soil thermodynamic to write restart files
1742    IF (hydrol_cwrr) THEN
1743       CALL thermosoil_finalize (kjit,    kjpindex, rest_id,   gtemp, &
1744            soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
1745    ELSE
1746       CALL thermosoilc_finalize (kjit,    kjpindex, rest_id,   gtemp, &
1747            soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
1748    END IF
1749
1750    !! 6. Add river routing to restart files 
1751    IF ( river_routing .AND. nbp_glo .GT. 1) THEN
1752       !! 6.1 Call river routing to write restart files
1753       CALL routing_finalize( kjit, kjpindex, rest_id, flood_frac, flood_res, streamfl_frac, fastr )
1754    ELSE
1755       !! 6.2 No routing, set variables to zero
1756       reinfiltration(:,:) = zero
1757       returnflow(:,:) = zero
1758       irrigation(:,:) = zero
1759           sed_deposition(:,:) = zero
1760           poc_deposition(:,:) = zero
1761           sed_deposition_d(:,:) = zero
1762           poc_deposition_d(:,:) = zero
1763       flood_frac(:) = zero
1764       streamfl_frac(:) = zero
1765       flood_res(:) = zero
1766       fastr(:) = zero   
1767    ENDIF
1768   
1769    !! 7. Call slowproc_main to add 'daily' and annual variables to restart file
1770    CALL slowproc_finalize (kjit,       kjpindex,  rest_id,  index,      &
1771                            njsc,       lai,       height,   veget,      &
1772                            frac_nobio, veget_max, reinf_slope,          &
1773                            co2_to_bm,  assim_param, frac_age, soiltile, &
1774                                                        litter_above, litter_below, carbon, DOC,     &
1775                                                        lignin_struc_above,lignin_struc_below,depth_deepsoil)
1776   
1777    IF (printlev_loc>=3) WRITE (numout,*) 'sechiba_finalize done'
1778   
1779  END SUBROUTINE sechiba_finalize
1780
1781 
1782!! ==============================================================================================================================\n
1783!! SUBROUTINE   : sechiba_init
1784!!
1785!>\BRIEF        Dynamic allocation of the variables, the dimensions of the
1786!! variables are determined by user-specified settings.
1787!!
1788!! DESCRIPTION  : The domain size (:: kjpindex) is used to allocate the correct
1789!! dimensions to all variables in sechiba. Depending on the variable, its
1790!! dimensions are also determined by the number of PFT's (::nvm), number of
1791!! soil types (::nstm), number of non-vegetative surface types (::nnobio),
1792!! number of soil levels (::ngrnd), number of soil layers in the hydrological
1793!! model (i.e. cwrr) (::nslm). Values for these variables are set in
1794!! constantes_soil.f90 and constantes_veg.f90.\n
1795!!
1796!! Memory is allocated for all Sechiba variables and new indexing tables
1797!! are build making use of both (::kjpij) and (::kjpindex). New indexing tables
1798!! are needed because a single pixel can contain several PFTs, soil types, etc.
1799!! The new indexing tables have separate indices for the different
1800!! PFTs, soil types, etc.\n
1801!!
1802!! RECENT CHANGE(S): None
1803!!
1804!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
1805!! variables. However, the routine allocates memory and builds new indexing
1806!! variables for later use.
1807!!
1808!! REFERENCE(S) : None
1809!!
1810!! FLOWCHART    : None
1811!! \n
1812!_ ================================================================================================================================
1813
1814  SUBROUTINE sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)
1815
1816!! 0.1 Input variables
1817 
1818    INTEGER(i_std), INTENT (in)                         :: kjit               !! Time step number (unitless)
1819    INTEGER(i_std), INTENT (in)                         :: kjpij              !! Total size of the un-compressed grid (unitless)
1820    INTEGER(i_std), INTENT (in)                         :: kjpindex           !! Domain size - terrestrial pixels only (unitless)
1821    INTEGER(i_std), INTENT (in)                         :: rest_id            !! _Restart_ file identifier (unitless)
1822    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: index              !! Indeces of the points on the map (unitless)
1823    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)     :: lalo               !! Geographical coordinates (latitude,longitude)
1824                                                                              !! for pixels (degrees)
1825!! 0.2 Output variables
1826
1827!! 0.3 Modified variables
1828
1829!! 0.4 Local variables
1830
1831    INTEGER(i_std)                                      :: ier                !! Check errors in memory allocation (unitless)
1832    INTEGER(i_std)                                      :: ji, jv             !! Indeces (unitless)
1833!_ ==============================================================================================================================
1834
1835!! 1. Initialize variables
1836   
1837    ! Dynamic allocation with user-specified dimensions on first call
1838    IF (l_first_sechiba) THEN
1839       l_first_sechiba=.FALSE.
1840    ELSE
1841       CALL ipslerr_p(3,'sechiba_init',' l_first_sechiba false . we stop ','','')
1842    ENDIF
1843
1844    !! Initialize local printlev
1845    printlev_loc=get_printlev('sechiba')
1846   
1847
1848    !! 1.1 Initialize 3D vegetation indexation table
1849    ALLOCATE (indexveg(kjpindex*nvm),stat=ier)
1850    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexveg','','')
1851
1852    ALLOCATE (indexlai(kjpindex*(nlai+1)),stat=ier)
1853    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlai','','')
1854
1855    ALLOCATE (indexsoil(kjpindex*nstm),stat=ier)
1856    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsoil','','')
1857
1858    ALLOCATE (indexnobio(kjpindex*nnobio),stat=ier)
1859    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnobio','','')
1860
1861    ALLOCATE (indexgrnd(kjpindex*ngrnd),stat=ier)
1862    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexgrnd','','')
1863
1864    ALLOCATE (indexsnow(kjpindex*nsnow),stat=ier)
1865    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsnow','','')
1866
1867    ALLOCATE (indexlayer(kjpindex*nslm),stat=ier)
1868    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlayer','','')
1869
1870    ALLOCATE (indexnslm(kjpindex*nslm),stat=ier)
1871    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnslm','','')
1872       
1873    ALLOCATE (indexalb(kjpindex*2),stat=ier)
1874    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexalb','','')
1875
1876    !! 1.2  Initialize 1D array allocation with restartable value
1877    ALLOCATE (flood_res(kjpindex),stat=ier)
1878    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for flood_res','','')
1879    flood_res(:) = undef_sechiba
1880
1881    ALLOCATE (fastr(kjpindex),stat=ier)
1882    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fastr','','')
1883    fastr(:) = undef_sechiba
1884
1885    ALLOCATE (flood_frac(kjpindex),stat=ier)
1886    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for flood_frac','','')
1887    flood_frac(:) = undef_sechiba
1888
1889    ALLOCATE (stream_frac(kjpindex),stat=ier)
1890    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for stream_frac','','')
1891    stream_frac(:) = undef_sechiba
1892
1893    ALLOCATE (streamfl_frac(kjpindex),stat=ier)
1894    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for streamfl_frac','','')
1895    streamfl_frac(:) = undef_sechiba
1896
1897    ALLOCATE (snow(kjpindex),stat=ier)
1898    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow','','')
1899    snow(:) = undef_sechiba
1900
1901    ALLOCATE (snow_age(kjpindex),stat=ier)
1902    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_age','','')
1903    snow_age(:) = undef_sechiba
1904
1905    ALLOCATE (drysoil_frac(kjpindex),stat=ier)
1906    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drysoil_frac','','')
1907
1908    ALLOCATE (rsol(kjpindex),stat=ier)
1909    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rsol','','')
1910
1911    ALLOCATE (evap_bare_lim(kjpindex),stat=ier)
1912    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim','','')
1913
1914    ALLOCATE (evapot(kjpindex),stat=ier)
1915    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot','','')
1916    evapot(:) = undef_sechiba
1917
1918    ALLOCATE (evapot_corr(kjpindex),stat=ier)
1919    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot_corr','','')
1920
1921    ALLOCATE (humrel(kjpindex,nvm),stat=ier)
1922    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for humrel','','')
1923    humrel(:,:) = undef_sechiba
1924
1925    ALLOCATE (vegstress(kjpindex,nvm),stat=ier)
1926    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vegstress','','')
1927    vegstress(:,:) = undef_sechiba
1928
1929    ALLOCATE (njsc(kjpindex),stat=ier)
1930    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for njsc','','')
1931    njsc(:)= undef_int
1932
1933    ALLOCATE (soiltile(kjpindex,nstm),stat=ier)
1934    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soiltile','','')
1935
1936    ALLOCATE (fraclut(kjpindex,nlut),stat=ier)
1937    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fraclut','','')
1938
1939    ALLOCATE (nwdFraclut(kjpindex,nlut),stat=ier)
1940    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for nwdFraclut','','')
1941
1942    ALLOCATE (reinf_slope(kjpindex),stat=ier)
1943    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope','','')
1944
1945    ALLOCATE (vbeta1(kjpindex),stat=ier)
1946    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta1','','')
1947
1948    ALLOCATE (vbeta4(kjpindex),stat=ier)
1949    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta4','','')
1950
1951    ALLOCATE (vbeta5(kjpindex),stat=ier)
1952    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta5','','')
1953
1954    ALLOCATE (soilcap(kjpindex),stat=ier)
1955    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilcap','','')
1956
1957    ALLOCATE (soilflx(kjpindex),stat=ier)
1958    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilflx','','')
1959
1960    ALLOCATE (temp_sol(kjpindex),stat=ier)
1961    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol','','')
1962    temp_sol(:) = undef_sechiba
1963
1964    ALLOCATE (qsurf(kjpindex),stat=ier)
1965    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsurf','','')
1966    qsurf(:) = undef_sechiba
1967
1968    !! 1.3 Initialize 2D array allocation with restartable value
1969    ALLOCATE (qsintveg(kjpindex,nvm),stat=ier)
1970    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintveg','','')
1971    qsintveg(:,:) = undef_sechiba
1972
1973    ALLOCATE (vbeta2(kjpindex,nvm),stat=ier)
1974    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta2','','')
1975
1976    ALLOCATE (vbeta3(kjpindex,nvm),stat=ier)
1977    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3','','')
1978
1979    ALLOCATE (vbeta3pot(kjpindex,nvm),stat=ier)
1980    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3pot','','')
1981
1982    ALLOCATE (gsmean(kjpindex,nvm),stat=ier)
1983    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gsmean','','')
1984
1985    ALLOCATE (cimean(kjpindex,nvm),stat=ier)
1986    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cimean','','')
1987
1988    ALLOCATE (gpp(kjpindex,nvm),stat=ier)
1989    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gpp','','')
1990    gpp(:,:) = undef_sechiba
1991
1992 
1993    ALLOCATE (temp_growth(kjpindex),stat=ier) 
1994    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_growth','','')
1995    temp_growth(:) = undef_sechiba 
1996
1997    ALLOCATE (veget(kjpindex,nvm),stat=ier)
1998    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget','','')
1999    veget(:,:)=undef_sechiba
2000
2001    ALLOCATE (veget_max(kjpindex,nvm),stat=ier)
2002    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget_max','','')
2003
2004    ALLOCATE (tot_bare_soil(kjpindex),stat=ier)
2005    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_bare_soil','','')
2006
2007    ALLOCATE (lai(kjpindex,nvm),stat=ier)
2008    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lai','','')
2009    lai(:,:)=undef_sechiba
2010
2011    ALLOCATE (frac_age(kjpindex,nvm,nleafages),stat=ier)
2012    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_age','','')
2013    frac_age(:,:,:)=undef_sechiba
2014
2015    ALLOCATE (height(kjpindex,nvm),stat=ier)
2016    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for height','','')
2017    height(:,:)=undef_sechiba
2018
2019    ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
2020    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_nobio','','')
2021    frac_nobio(:,:) = undef_sechiba
2022
2023    ALLOCATE (snow_nobio(kjpindex,nnobio),stat=ier)
2024    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio','','')
2025    snow_nobio(:,:) = undef_sechiba
2026
2027    ALLOCATE (snow_nobio_age(kjpindex,nnobio),stat=ier)
2028    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio_age','','')
2029    snow_nobio_age(:,:) = undef_sechiba
2030
2031    ALLOCATE (assim_param(kjpindex,nvm,npco2),stat=ier)
2032    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for assim_param','','')
2033
2034    !! 1.4 Initialize 1D array allocation
2035    ALLOCATE (vevapflo(kjpindex),stat=ier)
2036    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapflo','','')
2037    vevapflo(:)=zero
2038
2039    ALLOCATE (vevapsno(kjpindex),stat=ier)
2040    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapsno','','')
2041
2042    ALLOCATE (vevapnu(kjpindex),stat=ier)
2043    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapnu','','')
2044
2045    ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
2046    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for totfrac_nobio','','')
2047
2048    ALLOCATE (floodout(kjpindex),stat=ier)
2049    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for floodout','','')
2050
2051    ALLOCATE (runoff(kjpindex),stat=ier)
2052    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for runoff','','')
2053
2054    ALLOCATE (wat_flux0(kjpindex,nstm),stat=ier)
2055    IF (ier.NE.0) THEN
2056       WRITE (numout,*) ' error in wat_flux0 allocation. We stop. We need kjpindex*nstm words = ',kjpindex*nstm
2057       STOP 'sechiba_init'
2058    END IF
2059 
2060    ALLOCATE (wat_flux(kjpindex,nslm, nstm),stat=ier)
2061    IF (ier.NE.0) THEN
2062       WRITE (numout,*) ' error in wat_flux allocation. We stop. We need kjpindex*nslm*nstm words = ',kjpindex*nslm*nstm
2063       STOP 'sechiba_init'
2064    END IF
2065
2066    ALLOCATE (drainage_per_soil(kjpindex,nstm),stat=ier)
2067    IF (ier.NE.0) THEN
2068       WRITE (numout,*) ' error in drainage_per_soil allocation. We stop. We need kjpindex*nstm words = ',kjpindex*nstm
2069       STOP 'sechiba_init'
2070    END IF
2071
2072    ALLOCATE (runoff_per_soil(kjpindex,nstm),stat=ier)
2073    IF (ier.NE.0) THEN
2074       WRITE (numout,*) ' error in runoff_per_soil allocation. We stop. We need kjpindex*nstm words = ',kjpindex*nstm
2075       STOP 'sechiba_init'
2076    END IF
2077
2078    ALLOCATE (litter_mc(kjpindex,nstm),stat=ier)
2079    IF (ier.NE.0) THEN
2080       WRITE (numout,*) ' error in litter_mc allocation. We stop. We need kjpindex*nstm words = ',kjpindex*nstm
2081       STOP 'sechiba_init'
2082    END IF
2083
2084    ALLOCATE (soil_mc(kjpindex,nslm,nstm),stat=ier)
2085    IF (ier.NE.0) THEN
2086       WRITE (numout,*) ' error in soil_mc allocation. We stop. We need kjpindex*nslm*nstm words = ',kjpindex*nslm*nstm
2087       STOP 'sechiba_init'
2088    END IF
2089
2090    ALLOCATE (drainage(kjpindex),stat=ier)
2091    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drainage','','')
2092
2093    ALLOCATE (returnflow(kjpindex,nflow),stat=ier)
2094    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for returnflow','','')
2095    returnflow(:,:) = zero
2096   
2097    ALLOCATE (DOC_EXP_agg(kjpindex,nexp,nflow),stat=ier)
2098    IF (ier.NE.0) THEN
2099       WRITE (numout,*) ' error in DOC_EXP_agg allocation. We stop. We need kjpindex*nexp*nflow words = ',kjpindex*nexp*nflow
2100       STOP 'sechiba_init'
2101    END IF
2102    DOC_EXP_agg(:,:,:) = zero   
2103
2104        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Variables for erosion module zhc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2105        ALLOCATE (laydep(nslm),stat=ier)
2106    IF (ier.NE.0) THEN
2107       WRITE (numout,*) ' error in laydep allocation. We stop. We need kjpindex words = ',nslm
2108       STOP 'sechiba_init'
2109    END IF
2110    laydep(:) = zero
2111       
2112        ALLOCATE (std_totsed(kjpindex),stat=ier)
2113    IF (ier.NE.0) THEN
2114       WRITE (numout,*) ' error in std_totsed allocation. We stop. We need kjpindex words = ',kjpindex
2115       STOP 'sechiba_init'
2116    END IF
2117    std_totsed(:) = zero
2118       
2119        ALLOCATE (gridarea(kjpindex),stat=ier)
2120    IF (ier.NE.0) THEN
2121       WRITE (numout,*) ' error in gridarea allocation. We stop. We need kjpindex words = ',kjpindex
2122       STOP 'sechiba_init'
2123    END IF
2124    gridarea(:) = zero
2125       
2126        ALLOCATE (effgrid_perc(kjpindex),stat=ier)
2127    IF (ier.NE.0) THEN
2128       WRITE (numout,*) ' error in effgrid_perc allocation. We stop. We need kjpindex words = ',kjpindex
2129       STOP 'sechiba_init'
2130    END IF
2131    effgrid_perc(:) = zero
2132       
2133        ALLOCATE (bulkdens(kjpindex),stat=ier)
2134    IF (ier.NE.0) THEN
2135       WRITE (numout,*) ' error in bulkdens allocation. We stop. We need kjpindex words = ',kjpindex
2136       STOP 'sechiba_init'
2137    END IF
2138    bulkdens(:) = zero
2139       
2140        ALLOCATE (textfrac(kjpindex,nctext),stat=ier)
2141    IF (ier.NE.0) THEN
2142       WRITE (numout,*) ' error in textfrac allocation. We stop. We need kjpindex words = ',kjpindex*nctext
2143       STOP 'sechiba_init'
2144    END IF
2145    textfrac(:,:) = zero
2146       
2147        ALLOCATE (cf_veget(kjpindex,nvm),stat=ier)
2148    IF (ier.NE.0) THEN
2149       WRITE (numout,*) ' error in cf_veget allocation. We stop. We need kjpindex words = ',kjpindex*nvm
2150       STOP 'sechiba_init'
2151    END IF
2152    cf_veget(:,:) = zero
2153       
2154        ALLOCATE (kf_litter(kjpindex,nvm),stat=ier)
2155    IF (ier.NE.0) THEN
2156       WRITE (numout,*) ' error in kf_litter allocation. We stop. We need kjpindex words = ',kjpindex*nvm
2157       STOP 'sechiba_init'
2158    END IF
2159    kf_litter(:,:) = zero
2160       
2161        ALLOCATE (kf_root(kjpindex,nvm),stat=ier)
2162    IF (ier.NE.0) THEN
2163       WRITE (numout,*) ' error in kf_root allocation. We stop. We need kjpindex words = ',kjpindex*nvm
2164       STOP 'sechiba_init'
2165    END IF
2166    kf_root(:,:) = zero
2167
2168        ALLOCATE (totrunoff_d(kjpindex),stat=ier)
2169    IF (ier.NE.0) THEN
2170       WRITE (numout,*) ' error in totrunoff_d allocation. We stop. We need kjpindex words = ',kjpindex
2171       STOP 'sechiba_init'
2172    END IF
2173    totrunoff_d(:) = zero                                         
2174       
2175        ALLOCATE (peakrunoff_d(kjpindex),stat=ier)
2176    IF (ier.NE.0) THEN
2177       WRITE (numout,*) ' error in peakrunoff_d allocation. We stop. We need kjpindex words = ',kjpindex
2178       STOP 'sechiba_init'
2179    END IF
2180    peakrunoff_d(:) = zero
2181       
2182        ALLOCATE (totrundra_d(kjpindex),stat=ier)
2183    IF (ier.NE.0) THEN
2184       WRITE (numout,*) ' error in totrundra_d allocation. We stop. We need kjpindex words = ',kjpindex
2185       STOP 'sechiba_init'
2186    END IF
2187    totrundra_d(:) = zero
2188
2189        ALLOCATE (peakrundra_d(kjpindex),stat=ier)
2190    IF (ier.NE.0) THEN
2191       WRITE (numout,*) ' error in peakrundra_d allocation. We stop. We need kjpindex words = ',kjpindex
2192       STOP 'sechiba_init'
2193    END IF
2194    peakrundra_d(:) = zero
2195       
2196        ALLOCATE (erocap(kjpindex,nvm),stat=ier)
2197    IF (ier.NE.0) THEN
2198       WRITE (numout,*) ' error in erocap allocation. We stop. We need kjpindex words = ',kjpindex*nvm
2199       STOP 'sechiba_init'
2200    END IF
2201    erocap(:,:) = zero
2202       
2203        ALLOCATE (erodep(kjpindex,nvm),stat=ier)
2204    IF (ier.NE.0) THEN
2205       WRITE (numout,*) ' error in erodep allocation. We stop. We need kjpindex words = ',kjpindex*nvm
2206       STOP 'sechiba_init'
2207    erodep(:,:) = zero
2208    END IF
2209       
2210        ALLOCATE (eroc(kjpindex,ncarb,nvm),stat=ier)
2211    IF (ier.NE.0) THEN
2212       WRITE (numout,*) ' error in eroc allocation. We stop. We need kjpindex words = ',kjpindex*ncarb*nvm
2213       STOP 'sechiba_init'
2214    END IF
2215    eroc(:,:,:) = zero
2216       
2217        ALLOCATE (tot_erocap(kjpindex),stat=ier)
2218    IF (ier.NE.0) THEN
2219       WRITE (numout,*) ' error in tot_erocap allocation. We stop. We need kjpindex words = ',kjpindex
2220       STOP 'sechiba_init'
2221    END IF
2222    tot_erocap(:) = zero
2223       
2224        ALLOCATE (ave_erodep(kjpindex),stat=ier)
2225    IF (ier.NE.0) THEN
2226       WRITE (numout,*) ' error in ave_erodep allocation. We stop. We need kjpindex words = ',kjpindex
2227       STOP 'sechiba_init'
2228    ave_erodep(:) = zero
2229    END IF
2230       
2231        ALLOCATE (tot_eroc(kjpindex,ncarb),stat=ier)
2232    IF (ier.NE.0) THEN
2233       WRITE (numout,*) ' error in tot_eroc allocation. We stop. We need kjpindex words = ',kjpindex*ncarb
2234       STOP 'sechiba_init'
2235    END IF
2236    tot_eroc(:,:) = zero
2237       
2238        ALLOCATE (POC_EXP_agg(kjpindex,ncarb),stat=ier)
2239    IF (ier.NE.0) THEN
2240       WRITE (numout,*) ' error in POC_EXP_agg allocation. We stop. We need kjpindex words = ',kjpindex*ncarb
2241       STOP 'sechiba_init'
2242    END IF
2243    POC_EXP_agg(:,:) = zero
2244       
2245        ALLOCATE (DOC_ERO_agg(kjpindex,ncarb),stat=ier)
2246    IF (ier.NE.0) THEN
2247       WRITE (numout,*) ' error in DOC_ERO_agg allocation. We stop. We need kjpindex words = ',kjpindex*ncarb
2248       STOP 'sechiba_init'
2249    END IF
2250    DOC_ERO_agg(:,:) = zero
2251       
2252        ALLOCATE (SED_EXP_agg(kjpindex,nctext),stat=ier)
2253    IF (ier.NE.0) THEN
2254       WRITE (numout,*) ' error in SED_EXP_agg allocation. We stop. We need kjpindex words = ',kjpindex*nctext
2255       STOP 'sechiba_init'
2256    END IF
2257    SED_EXP_agg(:,:) = zero
2258       
2259        ALLOCATE (rootmass(kjpindex,nvm,nparts,nelements),stat=ier)
2260    IF (ier.NE.0) THEN
2261       WRITE (numout,*) ' error in rootmass allocation. We stop. We need kjpindex words = ',kjpindex*nvm*nparts*nelements
2262       STOP 'sechiba_init'
2263    END IF
2264    rootmass(:,:,:,:) = zero
2265       
2266        ALLOCATE (litter_above(kjpindex,nlitt,nvm,nelements),stat=ier)
2267    IF (ier.NE.0) THEN
2268       WRITE (numout,*) ' error in litter_above allocation. We stop. We need kjpindex words = ',kjpindex*nlitt*nvm*nelements
2269       STOP 'sechiba_init'
2270    END IF
2271    litter_above(:,:,:,:) = zero
2272       
2273        ALLOCATE (litter_below(kjpindex,nlitt,nvm,nslmd,nelements),stat=ier)
2274    IF (ier.NE.0) THEN
2275       WRITE (numout,*) ' error in litter_below allocation. We stop. We need kjpindex words = ',kjpindex*nlitt*nvm*(nslmd)*nelements
2276       STOP 'sechiba_init'
2277    END IF
2278    litter_below(:,:,:,:,:) = zero
2279       
2280        ALLOCATE (carbon(kjpindex,ncarb,nvm,nslmd),stat=ier)
2281    IF (ier.NE.0) THEN
2282       WRITE (numout,*) ' error in  allocation. We stop. We need kjpindex words = ',kjpindex*ncarb*nvm*(nslmd)
2283       STOP 'sechiba_init'
2284    END IF
2285    carbon(:,:,:,:) = zero
2286       
2287        ALLOCATE (DOC(kjpindex,nvm,nslmd,ndoc,npool,nelements),stat=ier)
2288    IF (ier.NE.0) THEN
2289       WRITE (numout,*) ' error in  allocation. We stop. We need kjpindex words = ',kjpindex*nvm*(nslmd)*ndoc*npool*nelements
2290       STOP 'sechiba_init'
2291    END IF
2292    DOC(:,:,:,:,:,:) = zero
2293       
2294        ALLOCATE(lignin_struc_above(kjpindex,nvm),stat=ier)
2295    IF (ier.NE.0) THEN
2296       WRITE(numout,*) 'Memory allocation error for lignin_struc_above. We stop. We need kjpindex*nvm words',kjpindex,nvm
2297       STOP 'sechiba_init'
2298    ENDIF
2299        lignin_struc_above(:,:) = zero
2300       
2301        ALLOCATE(lignin_struc_below(kjpindex,nvm,nslmd),stat=ier)
2302    IF (ier.NE.0) THEN
2303       WRITE(numout,*) 'Memory allocation error for lignin_struc_below. We stop. We need kjpindex*nvm*nlevs words',kjpindex,nvm,nslmd
2304       STOP 'sechiba_init'
2305    ENDIF
2306        lignin_struc_below(:,:,:) = zero
2307       
2308        ALLOCATE (depth_deepsoil(kjpindex,nvm),stat=ier)
2309    IF (ier.NE.0) THEN
2310       WRITE (numout,*) ' error in  allocation depth_deepsoil. We stop. We need kjpindex words = ',kjpindex*nvm
2311       STOP 'sechiba_init'
2312    END IF
2313    depth_deepsoil(:,:) = zero
2314       
2315        ALLOCATE (seddep_rate(kjpindex,nvm),stat=ier)
2316    IF (ier.NE.0) THEN
2317       WRITE (numout,*) ' error in  allocation seddep_rate. We stop. We need kjpindex words = ',kjpindex*nvm
2318       STOP 'sechiba_init'
2319    END IF
2320    seddep_rate(:,:) = zero
2321        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Variables for erosion module zhc !!!!!!!!!!!!!!!!!!!!!!
2322       
2323    ALLOCATE (reinfiltration(kjpindex,nflow),stat=ier)
2324    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinfiltration','','')
2325    reinfiltration(:,:) = zero
2326
2327    ALLOCATE (irrigation(kjpindex,nflow),stat=ier)
2328    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for irrigation','','')
2329    irrigation(:,:) = zero
2330       
2331        ALLOCATE (sed_deposition(kjpindex,nctext),stat=ier)
2332    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for sed_deposition','','')
2333    sed_deposition(:,:) = zero
2334       
2335        ALLOCATE (poc_deposition(kjpindex,ncarb),stat=ier)
2336    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for poc_deposition','','')
2337    poc_deposition(:,:) = zero
2338       
2339        ALLOCATE (sed_deposition_d(kjpindex,nctext),stat=ier)
2340    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for sed_deposition_d','','')
2341    sed_deposition_d(:,:) = zero
2342       
2343        ALLOCATE (poc_deposition_d(kjpindex,ncarb),stat=ier)
2344    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for poc_deposition_d','','')
2345    poc_deposition_d(:,:) = zero
2346   
2347    ALLOCATE (DOC_to_topsoil(kjpindex,nflow),stat=ier)
2348    IF (ier.NE.0) THEN
2349       WRITE (numout,*) ' error in DOC_to_topsoil allocation. We stop. We need kjpindex*nflow words = ',kjpindex*nflow
2350       STOP 'sechiba_init'
2351    END IF
2352    DOC_to_topsoil(:,:) = zero   
2353   
2354    ALLOCATE (DOC_to_subsoil(kjpindex,nflow),stat=ier)
2355    IF (ier.NE.0) THEN
2356       WRITE (numout,*) ' error in DOC_to_subsoil allocation. We stop. We need kjpindex*nflow words = ',kjpindex*nflow
2357       STOP 'sechiba_init'
2358    END IF
2359    DOC_to_subsoil(:,:) = zero
2360
2361    ALLOCATE (precip2canopy(kjpindex,nvm),stat=ier)
2362    IF (ier.NE.0) THEN
2363       WRITE (numout,*) ' error in precip2canopy allocation. We stop. We need kjpindex*nvm words = ',kjpindex*nvm
2364       STOP 'sechiba_init'
2365    END IF
2366
2367    ALLOCATE (precip2ground(kjpindex,nvm),stat=ier)
2368    IF (ier.NE.0) THEN
2369       WRITE (numout,*) ' error in precip2ground allocation. We stop. We need kjpindex*nvm words = ',kjpindex*nvm
2370       STOP 'sechiba_init'
2371    END IF
2372
2373    ALLOCATE (canopy2ground(kjpindex,nvm),stat=ier)
2374    IF (ier.NE.0) THEN
2375       WRITE (numout,*) ' error in canopy2ground allocation. We stop. We need kjpindex*nvm words = ',kjpindex*nvm
2376       STOP 'sechiba_init'
2377    END IF
2378   
2379    ALLOCATE (z0h(kjpindex),stat=ier)
2380    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0h','','')
2381
2382    ALLOCATE (z0m(kjpindex),stat=ier)
2383    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0m','','')
2384
2385    ALLOCATE (roughheight(kjpindex),stat=ier)
2386    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for roughheight','','')
2387
2388    ALLOCATE (emis(kjpindex),stat=ier)
2389    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for emis','','')
2390
2391    ALLOCATE (tot_melt(kjpindex),stat=ier)
2392    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_melt','','')
2393
2394    ALLOCATE (vbeta(kjpindex),stat=ier)
2395    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta','','')
2396
2397    ALLOCATE (fusion(kjpindex),stat=ier)
2398    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fusion','','')
2399
2400    ALLOCATE (rau(kjpindex),stat=ier)
2401    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rau','','')
2402
2403    ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
2404    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for deadleaf_cover','','')
2405
2406    ALLOCATE (stempdiag(kjpindex, nslm),stat=ier)
2407    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for stempdiag','','')
2408
2409    ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
2410    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for co2_flux','','')
2411    co2_flux(:,:)=zero
2412
2413    ALLOCATE (co2_to_bm(kjpindex,nvm),stat=ier)
2414    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for co2_to_bm','','')
2415
2416    ALLOCATE (shumdiag(kjpindex,nslm),stat=ier)
2417    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag','','')
2418   
2419    ALLOCATE (shumdiag_perma(kjpindex,nslm),stat=ier)
2420    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag_perma','','')
2421
2422    ALLOCATE (litterhumdiag(kjpindex),stat=ier)
2423    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for litterhumdiag','','')
2424
2425    ALLOCATE (ptnlev1(kjpindex),stat=ier)
2426    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ptnlev1','','')
2427
2428    ALLOCATE (k_litt(kjpindex),stat=ier)
2429    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for k_litt','','')
2430
2431    !! 1.5 Initialize 2D array allocation
2432    ALLOCATE (vevapwet(kjpindex,nvm),stat=ier)
2433    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapwet','','')
2434    vevapwet(:,:)=undef_sechiba
2435
2436    ALLOCATE (transpir(kjpindex,nvm),stat=ier)
2437    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpir','','')
2438
2439    ALLOCATE (transpot(kjpindex,nvm),stat=ier)
2440    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpot','','')
2441
2442    ALLOCATE (qsintmax(kjpindex,nvm),stat=ier)
2443    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintmax','','')
2444
2445    ALLOCATE (rveget(kjpindex,nvm),stat=ier)
2446    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rveget','','')
2447
2448    ALLOCATE (rstruct(kjpindex,nvm),stat=ier)
2449    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rstruct','','')
2450
2451    ALLOCATE (pgflux(kjpindex),stat=ier)
2452    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for pgflux','','')
2453    pgflux(:)= 0.0
2454
2455    ALLOCATE (cgrnd_snow(kjpindex,nsnow),stat=ier)
2456    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cgrnd_snow','','')
2457    cgrnd_snow(:,:) = 0
2458
2459    ALLOCATE (dgrnd_snow(kjpindex,nsnow),stat=ier)
2460    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for dgrnd_snow','','')
2461    dgrnd_snow(:,:) = 0
2462
2463    ALLOCATE (lambda_snow(kjpindex),stat=ier)
2464    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lambda_snow','','')
2465    lambda_snow(:) = 0
2466
2467    ALLOCATE (temp_sol_add(kjpindex),stat=ier)
2468    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol_add','','')
2469
2470    ALLOCATE (gtemp(kjpindex),stat=ier)
2471    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gtemp','','')
2472
2473    ALLOCATE (frac_snow_veg(kjpindex),stat=ier)
2474    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_snow_veg','','')
2475
2476    ALLOCATE (frac_snow_nobio(kjpindex,nnobio),stat=ier)
2477    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_snow_nobio','','')
2478
2479    ALLOCATE (snowrho(kjpindex,nsnow),stat=ier)
2480    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowrho','','')
2481
2482    ALLOCATE (snowheat(kjpindex,nsnow),stat=ier)
2483    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowheat','','')
2484
2485    ALLOCATE (snowgrain(kjpindex,nsnow),stat=ier)
2486    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowgrain','','')
2487
2488    ALLOCATE (snowtemp(kjpindex,nsnow),stat=ier)
2489    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowtemp','','')
2490
2491    ALLOCATE (snowdz(kjpindex,nsnow),stat=ier)
2492    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowdz','','')
2493
2494    ALLOCATE (mc_layh(kjpindex, nslm),stat=ier)
2495    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mc_layh','','')
2496
2497    ALLOCATE (mcl_layh(kjpindex, nslm),stat=ier)
2498    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcl_layh','','')
2499
2500    ALLOCATE (soilmoist(kjpindex, nslm),stat=ier)
2501    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilmoist','','')
2502
2503    !! 1.6 Initialize indexing table for the vegetation fields.
2504    ! In SECHIBA we work on reduced grids but to store in the full 3D filed vegetation variable
2505    ! we need another index table : indexveg, indexsoil, indexnobio and indexgrnd
2506    DO ji = 1, kjpindex
2507       !
2508       DO jv = 1, nlai+1
2509          indexlai((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2510       ENDDO
2511       !
2512       DO jv = 1, nvm
2513          indexveg((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2514       ENDDO
2515       !     
2516       DO jv = 1, nstm
2517          indexsoil((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2518       ENDDO
2519       !     
2520       DO jv = 1, nnobio
2521          indexnobio((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2522       ENDDO
2523       !
2524       DO jv = 1, ngrnd
2525          indexgrnd((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2526       ENDDO
2527       !
2528       DO jv = 1, nsnow
2529          indexsnow((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
2530       ENDDO
2531
2532       DO jv = 1, nslm
2533          indexnslm((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
2534       ENDDO
2535
2536       DO jv = 1, nslm
2537          indexlayer((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2538       ENDDO
2539       !
2540       DO jv = 1, 2
2541          indexalb((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2542       ENDDO
2543       !
2544    ENDDO
2545
2546!! 2. Read the default value that will be put into variable which are not in the restart file
2547    CALL ioget_expval(val_exp)
2548   
2549    IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_init done '
2550
2551  END SUBROUTINE sechiba_init
2552 
2553
2554!! ==============================================================================================================================\n
2555!! SUBROUTINE   : sechiba_clear
2556!!
2557!>\BRIEF        Deallocate memory of sechiba's variables
2558!!
2559!! DESCRIPTION  : None
2560!!
2561!! RECENT CHANGE(S): None
2562!!
2563!! MAIN OUTPUT VARIABLE(S): None
2564!!
2565!! REFERENCE(S) : None
2566!!
2567!! FLOWCHART    : None
2568!! \n
2569!_ ================================================================================================================================
2570
2571  SUBROUTINE sechiba_clear()
2572
2573!! 1. Initialize first run
2574    l_first_sechiba=.TRUE.
2575
2576!! 2. Deallocate dynamic variables of sechiba
2577    IF ( ALLOCATED (indexveg)) DEALLOCATE (indexveg)
2578    IF ( ALLOCATED (indexlai)) DEALLOCATE (indexlai)
2579    IF ( ALLOCATED (indexsoil)) DEALLOCATE (indexsoil)
2580    IF ( ALLOCATED (indexnobio)) DEALLOCATE (indexnobio)
2581    IF ( ALLOCATED (indexsnow)) DEALLOCATE (indexsnow)
2582    IF ( ALLOCATED (indexgrnd)) DEALLOCATE (indexgrnd)
2583    IF ( ALLOCATED (indexlayer)) DEALLOCATE (indexlayer)
2584    IF ( ALLOCATED (indexnslm)) DEALLOCATE (indexnslm)                                                                                             
2585    IF ( ALLOCATED (indexalb)) DEALLOCATE (indexalb)
2586    IF ( ALLOCATED (flood_res)) DEALLOCATE (flood_res)
2587    IF ( ALLOCATED (fastr)) DEALLOCATE (fastr)
2588    IF ( ALLOCATED (flood_frac)) DEALLOCATE (flood_frac)
2589    IF ( ALLOCATED (streamfl_frac)) DEALLOCATE (streamfl_frac)
2590    IF ( ALLOCATED (snow)) DEALLOCATE (snow)
2591    IF ( ALLOCATED (snow_age)) DEALLOCATE (snow_age)
2592    IF ( ALLOCATED (drysoil_frac)) DEALLOCATE (drysoil_frac)
2593    IF ( ALLOCATED (rsol)) DEALLOCATE (rsol)
2594    IF ( ALLOCATED (evap_bare_lim)) DEALLOCATE (evap_bare_lim)
2595    IF ( ALLOCATED (evapot)) DEALLOCATE (evapot)
2596    IF ( ALLOCATED (evapot_corr)) DEALLOCATE (evapot_corr)
2597    IF ( ALLOCATED (humrel)) DEALLOCATE (humrel)
2598    IF ( ALLOCATED (vegstress)) DEALLOCATE (vegstress)
2599    IF ( ALLOCATED (soiltile)) DEALLOCATE (soiltile)
2600    IF ( ALLOCATED (fraclut)) DEALLOCATE (fraclut)
2601    IF ( ALLOCATED (nwdFraclut)) DEALLOCATE (nwdFraclut)
2602    IF ( ALLOCATED (njsc)) DEALLOCATE (njsc)
2603    IF ( ALLOCATED (reinf_slope)) DEALLOCATE (reinf_slope)
2604    IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1)
2605    IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4)
2606    IF ( ALLOCATED (vbeta5)) DEALLOCATE (vbeta5)
2607    IF ( ALLOCATED (soilcap)) DEALLOCATE (soilcap)
2608    IF ( ALLOCATED (soilflx)) DEALLOCATE (soilflx)
2609    IF ( ALLOCATED (temp_sol)) DEALLOCATE (temp_sol)
2610    IF ( ALLOCATED (qsurf)) DEALLOCATE (qsurf)
2611    IF ( ALLOCATED (qsintveg)) DEALLOCATE (qsintveg)
2612    IF ( ALLOCATED (vbeta2))  DEALLOCATE (vbeta2)
2613    IF ( ALLOCATED (vbeta3)) DEALLOCATE (vbeta3)
2614    IF ( ALLOCATED (vbeta3pot)) DEALLOCATE (vbeta3pot)
2615    IF ( ALLOCATED (gsmean)) DEALLOCATE (gsmean)
2616    IF ( ALLOCATED (cimean)) DEALLOCATE (cimean)
2617    IF ( ALLOCATED (gpp)) DEALLOCATE (gpp)
2618    IF ( ALLOCATED (temp_growth)) DEALLOCATE (temp_growth) 
2619    IF ( ALLOCATED (veget)) DEALLOCATE (veget)
2620    IF ( ALLOCATED (veget_max)) DEALLOCATE (veget_max)
2621    IF ( ALLOCATED (tot_bare_soil)) DEALLOCATE (tot_bare_soil)
2622    IF ( ALLOCATED (lai)) DEALLOCATE (lai)
2623    IF ( ALLOCATED (frac_age)) DEALLOCATE (frac_age)
2624    IF ( ALLOCATED (height)) DEALLOCATE (height)
2625    IF ( ALLOCATED (roughheight)) DEALLOCATE (roughheight)
2626    IF ( ALLOCATED (frac_nobio)) DEALLOCATE (frac_nobio)
2627    IF ( ALLOCATED (snow_nobio)) DEALLOCATE (snow_nobio)
2628    IF ( ALLOCATED (snow_nobio_age)) DEALLOCATE (snow_nobio_age)
2629    IF ( ALLOCATED (assim_param)) DEALLOCATE (assim_param)
2630    IF ( ALLOCATED (vevapflo)) DEALLOCATE (vevapflo)
2631    IF ( ALLOCATED (vevapsno)) DEALLOCATE (vevapsno)
2632    IF ( ALLOCATED (vevapnu)) DEALLOCATE (vevapnu)
2633    IF ( ALLOCATED (totfrac_nobio)) DEALLOCATE (totfrac_nobio)
2634    IF ( ALLOCATED (floodout)) DEALLOCATE (floodout)
2635    IF ( ALLOCATED (runoff)) DEALLOCATE (runoff)
2636    IF ( ALLOCATED (drainage)) DEALLOCATE (drainage)
2637    IF ( ALLOCATED (reinfiltration)) DEALLOCATE (reinfiltration)
2638    IF ( ALLOCATED (DOC_EXP_agg)) DEALLOCATE (DOC_EXP_agg)
2639        !!!!!!!!!!!!!Variables for erosion ZHC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
2640    IF ( ALLOCATED (laydep)) DEALLOCATE (laydep)               
2641        IF ( ALLOCATED (bulkdens)) DEALLOCATE (bulkdens)
2642        IF ( ALLOCATED (textfrac)) DEALLOCATE (textfrac) 
2643        IF ( ALLOCATED (kf_litter)) DEALLOCATE (kf_litter)
2644        IF ( ALLOCATED (cf_veget)) DEALLOCATE (cf_veget)       
2645        IF ( ALLOCATED (kf_root)) DEALLOCATE (kf_root) 
2646        IF ( ALLOCATED (std_totsed)) DEALLOCATE (std_totsed)   
2647        IF ( ALLOCATED (gridarea)) DEALLOCATE (gridarea)       
2648        IF ( ALLOCATED (effgrid_perc)) DEALLOCATE (effgrid_perc)       
2649        IF ( ALLOCATED (totrunoff_d)) DEALLOCATE (totrunoff_d) 
2650        IF ( ALLOCATED (peakrunoff_d)) DEALLOCATE (peakrunoff_d)               
2651        IF ( ALLOCATED (totrundra_d)) DEALLOCATE (totrundra_d) 
2652        IF ( ALLOCATED (peakrundra_d)) DEALLOCATE (peakrundra_d)               
2653        IF ( ALLOCATED (erocap)) DEALLOCATE (erocap)   
2654        IF ( ALLOCATED (eroc)) DEALLOCATE (eroc)               
2655        IF ( ALLOCATED (erodep)) DEALLOCATE (erodep)   
2656        IF ( ALLOCATED (tot_erocap)) DEALLOCATE (tot_erocap)   
2657        IF ( ALLOCATED (tot_eroc)) DEALLOCATE (tot_eroc)       
2658        IF ( ALLOCATED (ave_erodep)) DEALLOCATE (ave_erodep)   
2659        IF ( ALLOCATED (POC_EXP_agg)) DEALLOCATE (POC_EXP_agg) 
2660        IF ( ALLOCATED (DOC_ERO_agg)) DEALLOCATE (DOC_ERO_agg)         
2661        IF ( ALLOCATED (SED_EXP_agg)) DEALLOCATE (SED_EXP_agg) 
2662        IF ( ALLOCATED (rootmass)) DEALLOCATE (rootmass)       
2663        IF ( ALLOCATED (litter_above)) DEALLOCATE (litter_above)       
2664        IF ( ALLOCATED (litter_below)) DEALLOCATE (litter_below)       
2665        IF ( ALLOCATED (carbon)) DEALLOCATE (carbon)   
2666        IF ( ALLOCATED (DOC)) DEALLOCATE (DOC) 
2667        IF ( ALLOCATED (lignin_struc_above)) DEALLOCATE (lignin_struc_above)   
2668        IF ( ALLOCATED (lignin_struc_below)) DEALLOCATE (lignin_struc_below)   
2669        IF ( ALLOCATED (depth_deepsoil)) DEALLOCATE (depth_deepsoil)   
2670        IF ( ALLOCATED (seddep_rate)) DEALLOCATE (seddep_rate)
2671!!!!!!!!!!!!!Variables for erosion ZHC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                                                                                                                           
2672    IF ( ALLOCATED (DOC_to_topsoil)) DEALLOCATE (DOC_to_topsoil) 
2673    IF ( ALLOCATED (DOC_to_subsoil)) DEALLOCATE (DOC_to_subsoil) 
2674    IF ( ALLOCATED (precip2canopy)) DEALLOCATE (precip2canopy)
2675    IF ( ALLOCATED (precip2ground)) DEALLOCATE (precip2ground)
2676    IF ( ALLOCATED (canopy2ground)) DEALLOCATE (canopy2ground)                                                                                                             
2677    IF ( ALLOCATED (irrigation)) DEALLOCATE (irrigation)
2678        IF ( ALLOCATED (sed_deposition)) DEALLOCATE (sed_deposition)
2679        IF ( ALLOCATED (poc_deposition)) DEALLOCATE (poc_deposition)
2680        IF ( ALLOCATED (sed_deposition_d)) DEALLOCATE (sed_deposition_d)
2681        IF ( ALLOCATED (poc_deposition_d)) DEALLOCATE (poc_deposition_d)
2682    IF ( ALLOCATED (tot_melt)) DEALLOCATE (tot_melt)
2683    IF ( ALLOCATED (vbeta)) DEALLOCATE (vbeta)
2684    IF ( ALLOCATED (fusion)) DEALLOCATE (fusion)
2685    IF ( ALLOCATED (rau)) DEALLOCATE (rau)
2686    IF ( ALLOCATED (deadleaf_cover)) DEALLOCATE (deadleaf_cover)
2687    IF ( ALLOCATED (stempdiag)) DEALLOCATE (stempdiag)
2688    IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
2689    IF ( ALLOCATED (co2_to_bm)) DEALLOCATE (co2_to_bm)
2690    IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
2691    IF ( ALLOCATED (shumdiag_perma)) DEALLOCATE (shumdiag_perma)
2692    IF ( ALLOCATED (litterhumdiag)) DEALLOCATE (litterhumdiag)
2693    IF ( ALLOCATED (ptnlev1)) DEALLOCATE (ptnlev1)
2694    IF ( ALLOCATED (k_litt)) DEALLOCATE (k_litt)
2695    IF ( ALLOCATED (vevapwet)) DEALLOCATE (vevapwet)
2696    IF ( ALLOCATED (transpir)) DEALLOCATE (transpir)
2697    IF ( ALLOCATED (transpot)) DEALLOCATE (transpot)
2698    IF ( ALLOCATED (qsintmax)) DEALLOCATE (qsintmax)
2699    IF ( ALLOCATED (rveget)) DEALLOCATE (rveget)
2700    IF ( ALLOCATED (rstruct)) DEALLOCATE (rstruct)
2701    IF ( ALLOCATED (frac_snow_veg)) DEALLOCATE (frac_snow_veg)
2702    IF ( ALLOCATED (frac_snow_nobio)) DEALLOCATE (frac_snow_nobio)
2703    IF ( ALLOCATED (snowrho)) DEALLOCATE (snowrho)
2704    IF ( ALLOCATED (snowgrain)) DEALLOCATE (snowgrain)
2705    IF ( ALLOCATED (snowtemp)) DEALLOCATE (snowtemp)
2706    IF ( ALLOCATED (snowdz)) DEALLOCATE (snowdz)
2707    IF ( ALLOCATED (snowheat)) DEALLOCATE (snowheat)
2708    IF ( ALLOCATED (cgrnd_snow)) DEALLOCATE (cgrnd_snow)
2709    IF ( ALLOCATED (dgrnd_snow)) DEALLOCATE (dgrnd_snow)
2710    IF ( ALLOCATED (lambda_snow)) DEALLOCATE(lambda_snow) 
2711    IF ( ALLOCATED (gtemp)) DEALLOCATE (gtemp)
2712    IF ( ALLOCATED (pgflux)) DEALLOCATE (pgflux)
2713    IF ( ALLOCATED (mc_layh)) DEALLOCATE (mc_layh)
2714    IF ( ALLOCATED (mcl_layh)) DEALLOCATE (mcl_layh)
2715    IF ( ALLOCATED (soilmoist)) DEALLOCATE (soilmoist)
2716
2717!! 3. Clear all allocated memory
2718    CALL pft_parameters_clear
2719    CALL slowproc_clear 
2720    CALL diffuco_clear 
2721    CALL enerbil_clear 
2722    IF ( hydrol_cwrr ) THEN
2723       CALL hydrol_clear 
2724       CALL thermosoil_clear
2725    ELSE
2726       CALL hydrolc_clear 
2727       CALL thermosoilc_clear
2728    ENDIF
2729    CALL condveg_clear 
2730    CALL routing_clear
2731
2732  END SUBROUTINE sechiba_clear
2733
2734
2735!! ==============================================================================================================================\n
2736!! SUBROUTINE   : sechiba_var_init
2737!!
2738!>\BRIEF        Calculate air density as a function of air temperature and
2739!! pressure for each terrestrial pixel.
2740!!
2741!! RECENT CHANGE(S): None
2742!!
2743!! MAIN OUTPUT VARIABLE(S): air density (::rau, kg m^{-3}).
2744!!
2745!! REFERENCE(S) : None
2746!!
2747!! FLOWCHART    : None
2748!! \n
2749!_ ================================================================================================================================
2750
2751  SUBROUTINE sechiba_var_init (kjpindex, rau, pb, temp_air) 
2752
2753!! 0.1 Input variables
2754
2755    INTEGER(i_std), INTENT (in)                    :: kjpindex        !! Domain size - terrestrial pixels only (unitless)
2756    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: pb              !! Surface pressure (hPa)
2757    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: temp_air        !! Air temperature (K)
2758   
2759!! 0.2 Output variables
2760
2761    REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rau             !! Air density @tex $(kg m^{-3})$ @endtex
2762
2763!! 0.3 Modified variables
2764
2765!! 0.4 Local variables
2766
2767    INTEGER(i_std)                                 :: ji              !! Indices (unitless)
2768!_ ================================================================================================================================
2769   
2770!! 1. Calculate intial air density (::rau)
2771   
2772    DO ji = 1,kjpindex
2773       rau(ji) = pa_par_hpa * pb(ji) / (cte_molr*temp_air(ji))
2774    END DO
2775
2776    IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_var_init done '
2777
2778  END SUBROUTINE sechiba_var_init
2779
2780
2781!! ==============================================================================================================================\n
2782!! SUBROUTINE   : sechiba_end
2783!!
2784!>\BRIEF        Swap old for newly calculated soil temperature.
2785!!
2786!! RECENT CHANGE(S): None
2787!!
2788!! MAIN OUTPUT VARIABLE(S): soil temperature (::temp_sol; K)
2789!!
2790!! REFERENCE(S) : None
2791!!
2792!! FLOWCHART    : None
2793!! \n
2794!! ================================================================================================================================
2795
2796  SUBROUTINE sechiba_end (kjpindex, temp_sol_new, temp_sol)
2797                         
2798
2799!! 0.1 Input variables
2800
2801    INTEGER(i_std), INTENT (in)                       :: kjpindex           !! Domain size - terrestrial pixels only (unitless)
2802    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: temp_sol_new       !! New soil temperature (K)
2803   
2804    !! 0.2 Output variables
2805    REAL(r_std),DIMENSION (kjpindex), INTENT (out)    :: temp_sol           !! Soil temperature (K)
2806
2807!_ ================================================================================================================================
2808   
2809!! 1. Swap temperature
2810
2811    temp_sol(:) = temp_sol_new(:)
2812   
2813    IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_end done '
2814
2815  END SUBROUTINE sechiba_end
2816
2817!! ==============================================================================================================================\n
2818!! SUBROUTINE   : sechiba_interface_orchidee_inca
2819!!
2820!>\BRIEF        make the interface between surface and atmospheric chemistry
2821!!
2822!! DESCRIPTION  : This subroutine is called from INCA, the atmospheric chemistry model. It is used to transfer variables from ORCHIDEE to INCA.
2823!!
2824!! RECENT CHANGE(S): move from chemistry module to be more generic (feb - 2017)
2825!!
2826!! MAIN OUTPUT VARIABLE(S): emission COV to be transport by orchidee to inca in fields_out array
2827!!
2828!! REFERENCE(S) : None
2829!!
2830!! FLOWCHART    : None
2831!! \n
2832!! ================================================================================================================================
2833  SUBROUTINE sechiba_interface_orchidee_inca( &
2834       nvm_out, veget_max_out, veget_frac_out, lai_out, snow_out, &
2835       field_out_names, fields_out, field_in_names, fields_in)
2836
2837
2838    INTEGER, INTENT(out)                      :: nvm_out            !! Number of vegetation types
2839    REAL(r_std), DIMENSION (:,:), INTENT(out) :: veget_max_out      !! Max. fraction of vegetation type (LAI -> infty)
2840    REAL(r_std), DIMENSION (:,:), INTENT(out) :: veget_frac_out     !! Fraction of vegetation type (unitless, 0-1) 
2841    REAL(r_std), DIMENSION (:,:), INTENT(out) :: lai_out            !! Surface foliere
2842    REAL(r_std), DIMENSION (:)  , INTENT(out) :: snow_out           !! Snow mass [Kg/m^2]
2843
2844    !
2845    ! Optional arguments
2846    !
2847    ! Names and fields for emission variables : to be transport by Orchidee to Inca
2848    CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_out_names
2849    REAL(r_std),DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: fields_out
2850    !
2851    ! Names and fields for deposit variables : to be transport from chemistry model by INCA to ORCHIDEE.
2852    CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_in_names
2853    REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(IN)    :: fields_in
2854
2855
2856    ! Variables always transmitted from sechiba to inca
2857    nvm_out = nvm 
2858    veget_max_out(:,:)  = veget_max(:,:) 
2859    veget_frac_out(:,:) = veget(:,:) 
2860    lai_out(:,:)  = lai(:,:) 
2861    snow_out(:)  = snow(:) 
2862
2863    ! Call chemistry_flux_interface if at least one of variables field_out_names or
2864    ! field_in_names is present in the argument list of sechiba_interface_orchidee_inca when called from inca.
2865    IF (PRESENT(field_out_names) .AND. .NOT. PRESENT(field_in_names)) THEN
2866       CALL chemistry_flux_interface(field_out_names=field_out_names, fields_out=fields_out)
2867    ELSE IF (.NOT. PRESENT(field_out_names) .AND. PRESENT(field_in_names)) THEN
2868       CALL chemistry_flux_interface(field_in_names=field_in_names, fields_in=fields_in)
2869    ELSE IF (PRESENT(field_out_names) .AND. PRESENT(field_in_names)) THEN
2870       CALL chemistry_flux_interface(field_out_names=field_out_names, fields_out=fields_out, &
2871            field_in_names=field_in_names, fields_in=fields_in)
2872    ENDIF
2873
2874  END SUBROUTINE sechiba_interface_orchidee_inca
2875
2876
2877END MODULE sechiba
2878
Note: See TracBrowser for help on using the repository browser.