Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_sechiba/sechiba.f90 @ 7852

Last change on this file since 7852 was 7615, checked in by sebastiaan.luyssaert, 2 years ago
Enhanced consistency of variable names: input has been changed in n_input throughout the code and the variable name vegstress introduced in sechiba is now also used in stomate. Enhnaced computational consistency: Pgap_cumul is used in stomate rather than recalculating it before calculating light_tran_to_floor_season. Edited getin_p while checking the code (but no real changes were made) and added several missing stomate and sechiba variables to age_class_distr to improve the 1+1=2 issue when comparing a model run with against a run without age classes. Finally: Write warning 10b in allocation to the history file rather than the out_orchidee file
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 166.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, hydrol_main,
11	!! 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): November 2020: It is possible to define soil hydraulic parameters from maps,
21	!! as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes).
22	!! Here, it leads to declare and allocate global variables.
23	!!
24	!! REFERENCE(S) : None
25	!!
26	!! SVN :
27	!! $HeadURL$
28	!! $Date$
29	!! $Revision$
30	!! \n
31	!_ ================================================================================================================================
32
33	MODULE sechiba
34
35	USE ioipsl
36	USE xios_orchidee
37
38	! modules used :
39	USE constantes
40	USE time, ONLY : one_day, dt_sechiba, FirstTsDay
41	USE solar, ONLY: solarang_noon
42	USE constantes_soil
43	USE pft_parameters
44	USE grid
45	USE structures
46	USE diffuco
47	USE condveg
48	USE enerbil
49	USE mleb
50	USE hydrol
51	USE hydraulic_arch
52	USE thermosoil
53	USE sechiba_io_p
54	USE slowproc
55	USE routing
56	USE ioipsl_para
57	USE chemistry
58	USE stomate_laieff
59	USE sapiens_lcchange, ONLY : check_veget
60
61	IMPLICIT NONE
62
63	PRIVATE
64	PUBLIC sechiba_main, sechiba_initialize, sechiba_clear, &
65	sechiba_interface_orchidee_inca, sechiba_xios_initialize
66
67	INTEGER(i_std), SAVE :: printlev_loc !! local printlev for this module
68	!$OMP THREADPRIVATE(printlev_loc)
69	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexveg !! indexing array for the 3D fields of vegetation
70	!$OMP THREADPRIVATE(indexveg)
71	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlai !! indexing array for the 3D fields of vegetation
72	!$OMP THREADPRIVATE(indexlai)
73	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnobio !! indexing array for the 3D fields of other surfaces (ice,
74	!! lakes, ...)
75	!$OMP THREADPRIVATE(indexnobio)
76	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsoil !! indexing array for the 3D fields of soil types (kjpindex*nstm)
77	!$OMP THREADPRIVATE(indexsoil)
78	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexgrnd !! indexing array for the 3D ground heat profiles (kjpindex*ngrnd)
79	!$OMP THREADPRIVATE(indexgrnd)
80	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlayer !! indexing array for the 3D fields of soil layers in CWRR (kjpindex*nslm)
81	!$OMP THREADPRIVATE(indexlayer)
82	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnslm !! indexing array for the 3D fields of diagnostic soil layers (kjpindex*nslm)
83	!$OMP THREADPRIVATE(indexnslm)
84	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexalb !! indexing array for the 2 fields of albedo
85	!$OMP THREADPRIVATE(indexalb)
86	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsnow !! indexing array for the 3D fields snow layers
87	!$OMP THREADPRIVATE(indexsnow)
88	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexcan !! indexing array for the level fields of the canopy
89	!$OMP THREADPRIVATE(indexcan)
90	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget !! Fraction of vegetation type (unitless, 0-1)
91	!$OMP THREADPRIVATE(veget)
92	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget_max !! Max. fraction of vegetation type (LAI -> infty, unitless)
93	!$OMP THREADPRIVATE(veget_max)
94	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_bare_soil !! Total evaporating bare soil fraction
95	!$OMP THREADPRIVATE(tot_bare_soil)
96	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: height !! Vegetation Height (m)
97	!$OMP THREADPRIVATE(height)
98	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: height_dom !! Dominant vegetation height (m)
99	!$OMP THREADPRIVATE(height_dom)
100	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
101	!! (unitless, 0-1)
102	!$OMP THREADPRIVATE(totfrac_nobio)
103	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: floodout !! Flow out of floodplains from hydrol
104	!$OMP THREADPRIVATE(floodout)
105	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: runoff !! Surface runoff calculated by hydrol
106	!! @tex $(kg m^{-2})$ @endtex
107	!$OMP THREADPRIVATE(runoff)
108	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: drainage !! Deep drainage calculatedd by hydrol
109	!! @tex $(kg m^{-2})$ @endtex
110	!$OMP THREADPRIVATE(drainage)
111	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: returnflow !! Water flow from lakes and swamps which returns to
112	!! the grid box @tex $(kg m^{-2})$ @endtex
113	!$OMP THREADPRIVATE(returnflow)
114	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: reinfiltration !! Routed water which returns into the soil
115	!$OMP THREADPRIVATE(reinfiltration)
116	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: irrigation !! Irrigation flux taken from the routing reservoirs and
117	!! being put into the upper layers of the soil
118	!! @tex $(kg m^{-2})$ @endtex
119	!$OMP THREADPRIVATE(irrigation)
120	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: emis !! Surface emissivity (unitless)
121	!$OMP THREADPRIVATE(emis)
122	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z0h !! Surface roughness for heat (m)
123	!$OMP THREADPRIVATE(z0h)
124	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z0m !! Surface roughness for momentum (m)
125	!$OMP THREADPRIVATE(z0m)
126	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: roughheight !! Effective height for roughness (m)
127	!$OMP THREADPRIVATE(roughheight)
128	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: reinf_slope !! slope coefficient (reinfiltration)
129	!$OMP THREADPRIVATE(reinf_slope)
130	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: ks !! Saturated soil conductivity (mm d^{-1})
131	!$OMP THREADPRIVATE(ks)
132	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: nvan !! Van Genushten n parameter (unitless)
133	!$OMP THREADPRIVATE(nvan)
134	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: avan !! Van Genushten alpha parameter (mm ^{-1})
135	!$OMP THREADPRIVATE(avan)
136	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcr !! Residual soil moisture (m^{3} m^{-3})
137	!$OMP THREADPRIVATE(mcr)
138	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcs !! Saturated soil moisture (m^{3} m^{-3})
139	!$OMP THREADPRIVATE(mcs)
140	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcfc !! Volumetric water content at field capacity (m^{3} m^{-3})
141	!$OMP THREADPRIVATE(mcfc)
142	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcw !! Volumetric water content at wilting point (m^{3} m^{-3})
143	!$OMP THREADPRIVATE(mcw)
144	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:,:):: us !! Water stress index for transpiration
145	!! (by soil layer and PFT) (0-1, unitless)
146	!$OMP THREADPRIVATE(us)
147	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shumdiag !! Mean relative soil moisture in the different levels used
148	!! by thermosoil.f90 (unitless, 0-1)
149	!$OMP THREADPRIVATE(shumdiag)
150	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shumdiag_perma !! Saturation degree of the soil
151	!$OMP THREADPRIVATE(shumdiag_perma)
152	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: k_litt !! litter cond.
153	!$OMP THREADPRIVATE(k_litt)
154	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: litterhumdiag !! Litter dryness factor (unitless, 0-1)
155	!$OMP THREADPRIVATE(litterhumdiag)
156	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: stempdiag !! Temperature which controls canopy evolution (K)
157	!$OMP THREADPRIVATE(stempdiag)
158	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: qsintveg !! Water on vegetation due to interception
159	!! @tex $(kg m^{-2})$ @endtex
160	!$OMP THREADPRIVATE(qsintveg)
161	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vbeta2 !! Interception resistance (unitless,0-1)
162	!$OMP THREADPRIVATE(vbeta2)
163	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vbeta3 !! Vegetation resistance (unitless,0-1)
164	!$OMP THREADPRIVATE(vbeta3)
165	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vbeta3pot !! Potential vegetation resistance
166	!$OMP THREADPRIVATE(vbeta3pot)
167	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gsmean !! Mean stomatal conductance for CO2 (mol m-2 s-1)
168	!$OMP THREADPRIVATE(gsmean)
169	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cimean !! STOMATE: mean intercellular CO2 concentration (ppm)
170	!$OMP THREADPRIVATE(cimean)
171	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vevapwet !! Interception loss over each PFT
172	!! @tex $(kg m^{-2} days^{-1})$ @endtex
173	!$OMP THREADPRIVATE(vevapwet)
174	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: transpir !! Transpiration @tex $(kg m^{-2} days^{-1})$ @endtex
175	!$OMP THREADPRIVATE(transpir)
176
177	! energy budget variables -----------------------------------------------
178	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: netrad !! Net radiation (W m^{-2})
179	!$OMP THREADPRIVATE(netrad)
180	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwabs !! LW radiation absorbed by the surface (W m^{-2})
181	!$OMP THREADPRIVATE(lwabs)
182	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwnet !! Net Long-wave radiation (W m^{-2})
183	!$OMP THREADPRIVATE(lwnet)
184	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: fluxsubli !! Energy of sublimation (mm day^{-1})
185	!$OMP THREADPRIVATE(fluxsubli)
186
187
188	! hydraulic stress variables -------------------------------------
189
190	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: transpir_supply !! Supply of water for transpiration @tex$$(mm dt^{-1})$ @endtex
191	!$OMP THREADPRIVATE(transpir_supply)
192	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vir_transpir_supply !! Supply of water for transpiration @tex$$(mm dt^{-1})$ @endtex
193	!$OMP THREADPRIVATE(vir_transpir_supply)
194	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: stressed !! Adjusted ecosystem functioning. Takes the unit of the variable
195	!! used as a proxy for waterstress
196	!$OMP THREADPRIVATE(stressed)
197	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: unstressed !! Initial ecosystem functioning after the first calculation and
198	!! before any recalculations. Takes the unit of the variable used
199	!! as a proxy for unstressed.
200	!$OMP THREADPRIVATE(unstressed)
201	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vessel_loss !! Conductivity lost due to cavitation in the xylem (no unit).
202	!$OMP THREADPRIVATE(vessel_loss)
203	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_frac !! Fraction of water transpired supplied by individual layers (no units)
204	!$OMP THREADPRIVATE(e_frac)
205
206	!$ REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vir_transpir_mod !! potential transpiration (transpot) divided by veget_max
207	!!$!$OMP THREADPRIVATE(vir_transpir_mod)
208	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: transpir_mod !! transpir divided by veget_max
209	!$OMP THREADPRIVATE(transpir_mod)
210
211	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:, :, :) :: transpir_column !! Supply of water for transpiration
212	!$OMP THREADPRIVATE(transpir_column)
213	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:, :, :) :: transpir_supply_column !! Supply of water for transpiration
214	!$OMP THREADPRIVATE(transpir_supply_column)
215	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:, :, :) :: transpir_mod_column !! Supply of water for transpiration
216	!$OMP THREADPRIVATE(transpir_mod_column)
217
218
219	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: transpot !! Potential transpiration (needed for irrigation)
220	!$OMP THREADPRIVATE(transpot)
221
222	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: qsintmax !! Maximum amount of water in the canopy interception
223	!! reservoir @tex $(kg m^{-2})$ @endtex
224	!$OMP THREADPRIVATE(qsintmax)
225	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: rveget !! Surface resistance for the vegetation
226	!! @tex $(s m^{-1})$ @endtex
227	!$OMP THREADPRIVATE(rveget)
228	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: rstruct !! Vegetation structural resistance
229	!$OMP THREADPRIVATE(rstruct)
230	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: warnings !! Holds a count of how many warnings we run into
231	!! of different types. It will get reset at the
232	!! end of each period since we don't wish to restart it.
233	!! Purely a technical diagnostic variable.
234	!$OMP THREADPRIVATE(warnings)
235	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: snow_nobio !! Snow mass of non-vegetative surfaces
236	!! @tex $(kg m^{-2})$ @endtex
237	!$OMP THREADPRIVATE(snow_nobio)
238	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: snow_nobio_age !! Snow age on non-vegetative surfaces (days)
239	!$OMP THREADPRIVATE(snow_nobio_age)
240	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: frac_nobio !! Fraction of non-vegetative surfaces (continental ice,
241	!! lakes, ...) (unitless, 0-1)
242	!$OMP THREADPRIVATE(frac_nobio)
243	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: assim_param !! vcmax, nue, and leaf N for photosynthesis for photosynthesis
244	!$OMP THREADPRIVATE(assim_param)
245	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: lai !! Surface foliaire
246	!$OMP THREADPRIVATE(lai)
247	TYPE(laieff_type), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: laieff_fit !! The parameters for fitting the effective
248	!! LAI function
249	!$OMP THREADPRIVATE(laieff_fit)
250	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gpp !! STOMATE: GPP. gC/m**2 of total area
251	!$OMP THREADPRIVATE(gpp)
252	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: temp_growth !! Growth temperature (C) - Is equal to t2m_month
253	!$OMP THREADPRIVATE(temp_growth)
254	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: humrel !! Relative soil moisture used in stomate to calculate plant water stress (0-1, unitless)
255	!$OMP THREADPRIVATE(humrel)
256	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vegstress !! Vegetation moisture stress (only for vegetation growth)
257	!$OMP THREADPRIVATE(vegstress)
258	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: frac_age !! Age efficacity from STOMATE for isoprene
259	!$OMP THREADPRIVATE(frac_age)
260	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: soiltile !! Fraction of each soil tile (0-1, unitless)
261	!$OMP THREADPRIVATE(soiltile)
262	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: fraclut !! Fraction of each landuse tile (0-1, unitless)
263	!$OMP THREADPRIVATE(fraclut)
264	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: nwdFraclut !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless)
265	!$OMP THREADPRIVATE(nwdFraclut)
266	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
267	!$OMP THREADPRIVATE(njsc)
268	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta1 !! Snow resistance
269	!$OMP THREADPRIVATE(vbeta1)
270	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta4 !! Bare soil resistance
271	!$OMP THREADPRIVATE(vbeta4)
272	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta5 !! Floodplains resistance
273	!$OMP THREADPRIVATE(vbeta5)
274	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: soilcap !!
275	!$OMP THREADPRIVATE(soilcap)
276	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: soilflx !!
277	!$OMP THREADPRIVATE(soilflx)
278	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol !! Soil temperature
279	!$OMP THREADPRIVATE(temp_sol)
280	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsurf !! near soil air moisture
281	!$OMP THREADPRIVATE(qsurf)
282	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: flood_res !! flood reservoir estimate
283	!$OMP THREADPRIVATE(flood_res)
284	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: flood_frac !! flooded fraction
285	!$OMP THREADPRIVATE(flood_frac)
286	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow !! Snow mass [Kg/m^2]
287	!$OMP THREADPRIVATE(snow)
288	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_age !! Snow age @tex ($d$) @endtex
289	!$OMP THREADPRIVATE(snow_age)
290	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: drysoil_frac !! Fraction of visibly (albedo) Dry soil (Between 0 and 1)
291	!$OMP THREADPRIVATE(drysoil_frac)
292	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: evap_bare_lim !! Bare soil stress
293	!$OMP THREADPRIVATE(evap_bare_lim)
294	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) ::swc !! Soil water content (a copy of mc) m3 m-3
295	!$OMP THREADPRIVATE(swc)
296	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) ::ksoil !! Soil conductivity (copy of k mm/d)
297	!$OMP THREADPRIVATE(ksoil)
298	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: evap_bare_lim_ns !! Bare soil stress
299	!$OMP THREADPRIVATE(evap_bare_lim_ns)
300	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: co2_to_bm !! virtual CO2 flux (gC/m**2 of average ground/s)
301	!$OMP THREADPRIVATE(co2_to_bm)
302	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: evapot !! Soil Potential Evaporation
303	!$OMP THREADPRIVATE(evapot)
304	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: evapot_corr !! Soil Potential Evaporation Correction (Milly 1992)
305	!$OMP THREADPRIVATE(evapot_corr)
306	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vevapflo !! Floodplains evaporation
307	!$OMP THREADPRIVATE(vevapflo)
308	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vevapsno !! Snow evaporation
309	!$OMP THREADPRIVATE(vevapsno)
310	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vevapnu !! Bare soil evaporation
311	!$OMP THREADPRIVATE(vevapnu)
312	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_melt !! Total melt
313	!$OMP THREADPRIVATE(tot_melt)
314	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta !! Resistance coefficient
315	!$OMP THREADPRIVATE(vbeta)
316	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: rau !! Density
317	!$OMP THREADPRIVATE(rau)
318	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: deadleaf_cover !! Fraction of soil covered by dead leaves
319	!$OMP THREADPRIVATE(deadleaf_cover)
320	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: ptnlev1 !! 1st level Different levels soil temperature
321	!$OMP THREADPRIVATE(ptnlev1)
322	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: laieff_isotrop !! Effective LAI
323	!$OMP THREADPRIVATE(laieff_isotrop)
324	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: Light_Abs_Tot !! Absorbed radiation per level for photosynthesis
325	!$OMP THREADPRIVATE(Light_Abs_Tot)
326	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: Light_Tran_Tot !! Transmitted radiation per level for photosynthesis
327	!$OMP THREADPRIVATE(Light_Tran_Tot)
328	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mc_layh !! Volumetric soil moisture for each layer in hydrol(liquid + ice) (m3/m3)
329	!$OMP THREADPRIVATE(mc_layh)
330	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mcl_layh !! Volumetric soil moisture for each layer in hydrol(liquid) (m3/m3)
331	!$OMP THREADPRIVATE(mcl_layh)
332	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: soilmoist !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
333	!$OMP THREADPRIVATE(soilmoist)
334	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: mc_layh_s !! Volumetric soil moisture for each layer in hydrol per soiltile (liquid + ice) (m3/m3)
335	!$OMP THREADPRIVATE(mc_layh_s)
336	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: mcl_layh_s !! Volumetric soil moisture for each layer in hydrol per soiltile (liquid) (m3/m3)
337	!$OMP THREADPRIVATE(mcl_layh_s)
338	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: soilmoist_s !! Total soil moisture content for each layer in hydrol per soiltile (liquid + ice) (mm)
339	!$OMP THREADPRIVATE(soilmoist_s)
340
341	! multi-layer variables --------------------------------
342	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: u_speed !! Canopy wind speed profile
343	!$OMP THREADPRIVATE(u_speed)
344	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: circ_class_biomass !! Stem diameter @tex $(m)$ @endtex
345	!$OMP THREADPRIVATE(circ_class_biomass)
346	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: circ_class_n !! Number of trees within each circumference
347	!$OMP THREADPRIVATE(circ_class_n)
348	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: loss_gain !! loss and gains in each PFT due to LCC
349	!$OMP THREADPRIVATE(loss_gain)
350	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget_max_new !! New year fraction of vegetation type (0-1, unitless)
351	!$OMP THREADPRIVATE(veget_max_new)
352	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: frac_nobio_new !! New year fraction of ice+lakes+cities+... (0-1, unitless)
353	!$OMP THREADPRIVATE(frac_nobio_new)
354	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: nlevels_loc !! The number of physical levels in the canopy
355	!$OMP THREADPRIVATE(nlevels_loc)
356
357	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: sum_veget_diff !! The is the difference between the total
358	!! vegetation fraction in the new and old
359	!! land cover maps. Useful for land cover
360	!! changes. [-]
361	!$OMP THREADPRIVATE(sum_veget_diff)
362
363	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:,:) :: z_array_out !! heights of tree levels from stomate
364	!$OMP THREADPRIVATE(z_array_out)
365
366	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: profile_vbeta3
367	!$OMP THREADPRIVATE(profile_vbeta3)
368
369	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: profile_rveget
370	!$OMP THREADPRIVATE(profile_rveget)
371
372	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: max_height_store !! Same as z_array, but one less dimension.
373	!! @tex $(m)$ @endte
374	!$OMP THREADPRIVATE(max_height_store)
375
376	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: delta_c13_assim !! C13 concentration in delta notation
377	!! @tex $ permille $ @endtex (per thousand)
378	!$OMP THREADPRIVATE(delta_c13_assim)
379
380	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: leaf_ci_out !! Ci Leaf internal CO2 concentration ()
381	!$OMP THREADPRIVATE(leaf_ci_out)
382
383	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gpp_day !! Number of time steps when there is gpp
384	!$OMP THREADPRIVATE(gpp_day)
385
386	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: lai_per_level !! The LAI per vertical level
387	!! @tex $(m^2 / m^2)$ @endtex
388	!$OMP THREADPRIVATE(lai_per_level)
389	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: frac_snow_pix !! The fraction of the whole pixel covered
390	!! by snow. This is computed from the above
391	!! two. @tex $-$ @endtex
392	!$OMP THREADPRIVATE(frac_snow_pix)
393
394
395	LOGICAL, SAVE :: l_first_sechiba = .TRUE. !! Flag controlling the intialisation (true/false)
396	!$OMP THREADPRIVATE(l_first_sechiba)
397
398	! Variables related to snow processes calculations
399	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: frac_snow_veg !! Snow cover fraction on vegetation (unitless)
400	!$OMP THREADPRIVATE(frac_snow_veg)
401	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: frac_snow_nobio !! Snow cover fraction on continental ice, lakes, etc (unitless)
402	!$OMP THREADPRIVATE(frac_snow_nobio)
403	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: snowrho !! snow density for each layer (Kg/m^3)
404	!$OMP THREADPRIVATE(snowrho)
405	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: snowheat !! snow heat content for each layer (J/m2)
406	!$OMP THREADPRIVATE(snowheat)
407	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: snowgrain !! snow grain size (m)
408	!$OMP THREADPRIVATE(snowgrain)
409	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: snowtemp !! snow temperature profile (K)
410	!$OMP THREADPRIVATE(snowtemp)
411	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: snowdz !! snow layer thickness (m)
412	!$OMP THREADPRIVATE(snowdz)
413	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: gtemp !! soil surface temperature
414	!$OMP THREADPRIVATE(gtemp)
415	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: pgflux !! net energy into snow pack
416	!$OMP THREADPRIVATE(pgflux)
417	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
418	!$OMP THREADPRIVATE(cgrnd_snow)
419	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
420	!$OMP THREADPRIVATE(dgrnd_snow)
421	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
422	!! from the first and second snow layers
423	!$OMP THREADPRIVATE(lambda_snow)
424	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
425	!$OMP THREADPRIVATE(temp_sol_add)
426	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: qsol_sat_new !! New saturated surface air moisture (kg kg^{-1})
427	!$OMP THREADPRIVATE(qsol_sat_new)
428	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: qair_new !! New specific humidity at lowest level (kg kg^{-1})
429	!$OMP THREADPRIVATE(qair_new)
430
431	!+++CHECK+++
432	! Variables defined in CN-CAN but no longer present in CN
433	LOGICAL, ALLOCATABLE, SAVE, DIMENSION(:) :: osfcmelt !! Indicate snow melting in each gridcell
434	!$OMP THREADPRIVATE(osfcmelt)
435	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: snowflx !! Snow flux (W m^{-2})
436	!$OMP THREADPRIVATE(snowflx)
437	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: snowcap !! Snow calorific capacity (J K^{-1])
438	!$OMP THREADPRIVATE(snowcap)
439	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: cgrnd_soil !! matrix coefficient for the computation of soil, from thermosoil
440	!$OMP THREADPRIVATE(cgrnd_soil)
441	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: dgrnd_soil !! matrix coefficient for the computation of soil, from thermosoil
442	!$OMP THREADPRIVATE(dgrnd_soil)
443	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: zdz1_soil !! numerical constant from thermosoil
444	!$OMP THREADPRIVATE(zdz1_soil)
445	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: zdz2_soil !! numerical constant from thermosoil
446	!$OMP THREADPRIVATE(zdz2_soil)
447	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: albedo_undersnow !! albedo under the snowpack
448	!$OMP THREADPRIVATE(albedo_undersnow)
449	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa_snow !! snow thermal conductivity
450	!$OMP THREADPRIVATE(pkappa_snow)
451	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: gthick !! soil surface layer thickness
452	!$OMP THREADPRIVATE(gthick)
453
454	!+++++++++++
455
456	! Variables related to carbon soil discretization
457	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tdeep !! Deep temperature profile (K)
458	!$OMP THREADPRIVATE(tdeep)
459	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hsdeep !! Deep soil humidity profile (unitless)
460	!$OMP THREADPRIVATE(hsdeep)
461	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: heat_Zimov !! Heating associated with decomposition [W/m**3 soil]
462	!$OMP THREADPRIVATE(heat_Zimov)
463	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: sfluxCH4_deep !! Surface flux of CH4 to atmosphere from permafrost
464	!$OMP THREADPRIVATE(sfluxCH4_deep)
465	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: sfluxCO2_deep !! Surface flux of CO2 to atmosphere from permafrost
466	!$OMP THREADPRIVATE(sfluxCO2_deep)
467	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: som_total !! total soil organic matter for use in thermal calcs (g/m**3)
468	!$OMP THREADPRIVATE(som_total)
469	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: altmax !! Maximul active layer thickness (m). Be careful, here active means non frozen.
470	!! Not related with the active soil carbon pool.
471	!$OMP THREADPRIVATE(altmax)
472	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: depth_organic_soil !! Depth at which there is still organic matter (m)
473	!$OMP THREADPRIVATE(depth_organic_soil)
474
475	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:,:) :: root_profile !! Normalized root length fraction in each soil layer
476	!! (0-1, unitless)
477	!$OMP THREADPRIVATE(root_profile)
478	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: root_depth !! Node and interface numbers at which the deepest roots
479	!! occur (1 to nslm, unitless)
480	!$OMP THREADPRIVATE(root_depth)
481	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coszang_noon !! Solar zenith angle at noon (dimensionless)
482	!$OMP THREADPRIVATE(coszang_noon)
483	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: Pgap_cumul !! The probability of finding a gap in the in canopy from the top
484	!! of the canopy to a given level (unitless, between 0-1)
485	!$OMP THREADPRIVATE(Pgap_cumul)
486
487	CONTAINS
488
489
490	!! =============================================================================================================================
491	!! SUBROUTINE: sechiba_xios_initialize
492	!!
493	!>\BRIEF Initialize xios dependant defintion before closing context defintion
494	!!
495	!! DESCRIPTION: Initialize xios dependant defintion before closing context defintion
496	!!
497	!! \n
498	!_ ==============================================================================================================================
499
500	SUBROUTINE sechiba_xios_initialize
501	LOGICAL :: lerr
502
503	IF (xios_orchidee_ok) THEN
504	lerr=xios_orchidee_setvar('min_sechiba',min_sechiba)
505	CALL slowproc_xios_initialize
506	CALL condveg_xios_initialize
507	CALL chemistry_xios_initialize
508	CALL thermosoil_xios_initialize
509	END IF
510	IF (printlev_loc>=3) WRITE(numout,*) 'End sechiba_xios_initialize'
511
512	END SUBROUTINE sechiba_xios_initialize
513
514
515
516
517	!! =============================================================================================================================
518	!! SUBROUTINE: sechiba_initialize
519	!!
520	!>\BRIEF Initialize all prinicipal modules by calling their "_initialize" subroutines
521	!!
522	!! DESCRIPTION: Initialize all prinicipal modules by calling their "_initialize" subroutines
523	!!
524	!! \n
525	!_ ==============================================================================================================================
526
527	SUBROUTINE sechiba_initialize( &
528	kjit, kjpij, kjpindex, index, &
529	lalo, contfrac, neighbours, resolution, zlev, &
530	u, v, qair, temp_air, &
531	petAcoef, peqAcoef, petBcoef, peqBcoef, &
532	precip_rain, precip_snow, lwdown, swnet, swdown, &
533	pb, rest_id, hist_id, hist2_id, &
534	rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
535	coastalflow, riverflow, tsol_rad, vevapp, qsurf_out, &
536	z0m_out, z0h_out, albedo, fluxsens, fluxlat, emis_out, &
537	temp_sol_new, tq_cdrag, coszang)
538
539	!! 0.1 Input variables
540	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
541	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
542	!! (unitless)
543	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
544	!! (unitless)
545	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
546	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (unitless)
547	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier (unitless)
548	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
549	!! (unitless)
550	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _History_ file identifier
551	!! (unitless)
552	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file
553	!! identifier (unitless)
554	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographic coordinates (latitude,longitude)
555	!! for grid cells (degrees)
556	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid
557	!! (unitless, 0-1)
558	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the pixels on the map.
559	!! Sechiba uses a reduced grid excluding oceans
560	!! ::index contains the indices of the
561	!! terrestrial pixels only! (unitless)
562	INTEGER(i_std), DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours !! Neighboring grid points if land!(unitless)
563	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size in x and y of the grid (m)
564
565	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed in direction u
566	!! @tex $(m.s^{-1})$ @endtex
567	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed in direction v
568	!! @tex $(m.s^{-1})$ @endtex
569	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
570	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
571	!! @tex $(kg kg^{-1})$ @endtex
572	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
573	!! @tex $(kg m^{-2})$ @endtex
574	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
575	!! @tex $(kg m^{-2})$ @endtex
576	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux
577	!! @tex $(W m^{-2})$ @endtex
578	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux
579	!! @tex $(W m^{-2})$ @endtex
580	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux
581	!! @tex $(W m^{-2})$ @endtex
582	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
583	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! Coefficients A for T from the Planetary
584	!! Boundary Layer
585	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! Coefficients A for q from the Planetary
586	!! Boundary Layer
587	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! Coefficients B for T from the Planetary
588	!! Boundary Layer
589	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! Coefficients B for q from the Planetary
590	!! Boundary Layer
591	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
592	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: coszang !! cosine of solar zenith angle
593
594
595	!! 0.2 Output variables
596	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: coastalflow !! Outflow on coastal points by small basins.
597	!! This is the water which flows in a disperse
598	!! way into the ocean
599	!! @tex $(kg dt_routing^{-1})$ @endtex
600	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: riverflow !! Outflow of the major rivers.
601	!! The flux will be located on the continental
602	!! grid but this should be a coastal point
603	!! @tex $(kg dt_routing^{-1})$ @endtex
604	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature
605	!! @tex $(W m^{-2})$ @endtex
606	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation
607	!! @tex $(kg m^{-2} days^{-1})$ @endtex
608
609	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf_out !! Surface specific humidity
610	!! @tex $(kg kg^{-1})$ @endtex
611	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m_out !! Surface roughness momentum (output diagnostic, m)
612	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h_out !! Surface roughness heat (output diagnostic, m)
613	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Surface albedo for visible and near-infrared
614	!! (unitless, 0-1)
615	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux
616	!! @tex $(W m^{-2})$ @endtex
617	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux
618	!! @tex $(W m^{-2})$ @endtex
619	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis_out !! Emissivity (unitless)
620	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New ground temperature (K)
621
622	!! 0.3 Modified
623	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: tq_cdrag !! Surface drag coefficient (-)
624
625	!! 0.4 Local variables
626	INTEGER(i_std) :: ji, jv, ilev !! Index (unitless)
627	REAL(r_std), DIMENSION(kjpindex) :: zmaxh_glo !! 2D field of constant soil depth (zmaxh) (m)
628	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
629	REAL(r_std),DIMENSION (kjpindex) :: epot_air !! Air potential energy (??J)
630	!+++CHECK+++
631	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
632	! when running a larger domain. Needs to be corrected when implementing
633	! a global use of the multi-layer energy budget.
634	!!$ REAL(r_std),DIMENSION(nlevels_tot) :: Light_Abs_Tot_mean!! total light absorption for a given canopy level
635	!!$ REAL(r_std),DIMENSION(nlevels_tot) :: Light_Alb_Tot_mean!! total albedo for a given level
636	!+++++++++++
637	INTEGER(i_std) :: init_config !! Identifer of the configuration used to
638	!! initialize stomate/or sechiba
639	REAL(r_std), DIMENSION(kjpindex,nslm,nvm) :: mc_layh_pft !! As mc_layh per pft
640	REAL(r_std), DIMENSION(kjpindex,nslm,nvm) :: mcl_layh_pft !! As mcl_layh per pft
641	REAL(r_std), DIMENSION(kjpindex,nslm,nvm) :: soilmoist_pft !! As soilmoist per pft
642
643	!_ ================================================================================================================================
644
645	!! Initialize local printlev
646	printlev_loc=get_printlev('sechiba')
647
648	IF (printlev_loc>=3) WRITE(numout,*) 'Start sechiba_initialize'
649
650	!! 1. Initialize variables on first call
651
652	!! 1.2 Initialize most of sechiba's variables
653	CALL sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)
654
655	!! 1.3 Initialize stomate's variables
656	CALL slowproc_initialize (kjit, kjpij, kjpindex, &
657	rest_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
658	index, indexveg, lalo, neighbours, &
659	resolution, contfrac, temp_air, coszang_noon, &
660	soiltile, reinf_slope, deadleaf_cover, assim_param, &
661	ks, nvan, avan, mcr, &
662	mcs, mcfc, mcw, &
663	frac_age, height, lai, veget, &
664	frac_nobio, njsc, veget_max, fraclut, &
665	nwdfraclut, tot_bare_soil,totfrac_nobio, qsintmax, &
666	temp_growth, circ_class_biomass, &
667	circ_class_n, lai_per_level,laieff_fit, &
668	z_array_out, max_height_store, &
669	som_total, heat_Zimov, altmax, depth_organic_soil,&
670	loss_gain, veget_max_new,frac_nobio_new, height_dom, &
671	Pgap_cumul)
672
673	!! 1.4 Initialize diffusion coefficients
674	CALL diffuco_initialize (kjit, kjpindex, index, &
675	rest_id, lalo, neighbours, resolution, &
676	rstruct, tq_cdrag)
677
678
679	!! 1.5 Initialize variables for the energy budget
680	IF (ok_mleb) THEN
681	!! Multi-layer energy budget
682	CALL mleb_initialize (kjit, kjpindex, rest_id, &
683	temp_air, qair, &
684	temp_sol, temp_sol_new, tsol_rad, &
685	evapot, evapot_corr, qsurf, fluxsens, &
686	fluxlat, vevapp, &
687	u_speed, z_array_out, &
688	max_height_store )
689	ELSE
690	!! Energy budget using enerbil module
691	CALL enerbil_initialize (kjit, kjpindex, index, rest_id, &
692	qair, &
693	temp_sol, temp_sol_new, tsol_rad, &
694	evapot, evapot_corr, qsurf, fluxsens, &
695	fluxlat, vevapp )
696	END IF
697
698	!! do we need to call enerbill here... it is called in enerbil_main, but we might needs these variables
699	!! for the multi-layer energy budget. Maybe we should put all this in the ok_mleb if statement? asla, MERGE
700	!! it crashes... investigate further
701	!!CALL enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis)
702
703
704
705	!! 1.7 Initialize remaining hydrological variables
706	CALL hydrol_initialize (ks, nvan, avan, mcr, mcs, mcfc, mcw, &
707	kjit, kjpindex, index, rest_id, &
708	njsc, soiltile, veget, veget_max, &
709	altmax, &
710	humrel, vegstress, drysoil_frac, &
711	shumdiag_perma, qsintveg, &
712	evap_bare_lim, evap_bare_lim_ns, snow, snow_age, snow_nobio, &
713	snow_nobio_age, snowrho, snowtemp, snowgrain, &
714	snowdz, snowheat, &
715	mc_layh, mcl_layh, soilmoist, mc_layh_s, mcl_layh_s, &
716	soilmoist_s, swc, ksoil, root_profile, &
717	us)
718
719	!! 1.9 Initialize surface parameters (emissivity, albedo and roughness)
720	!+++CHECK+++
721	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
722	! when running a larger domain. Needs to be corrected when implementing
723	! a global use of the multi-layer energy budget.
724	CALL condveg_initialize (kjit, kjpindex, index, rest_id, &
725	lalo, neighbours, resolution, contfrac, &
726	veget, veget_max, frac_nobio, totfrac_nobio, &
727	zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
728	drysoil_frac, height, height_dom, snowdz,snowrho, tot_bare_soil, &
729	temp_air, pb, u, v, &
730	lai, &
731	emis, albedo, z0m, z0h, roughheight, &
732	frac_snow_veg,frac_snow_nobio, &
733	coszang, &
734	Light_Abs_Tot, Light_Tran_Tot, laieff_fit, &
735	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
736	laieff_isotrop)
737	!+++++++++++
738
739	!! 1.10 Initialization of soil thermodynamics
740	DO jv = 1,nvm
741	mc_layh_pft(:,:,jv) = mc_layh_s(:,:,pref_soil_veg(jv))
742	mcl_layh_pft(:,:,jv) = mcl_layh_s(:,:,pref_soil_veg(jv))
743	soilmoist_pft(:,:,jv) = soilmoist_s(:,:,pref_soil_veg(jv))
744	END DO
745	CALL thermosoil_initialize (kjit, kjpindex, rest_id, &
746	temp_sol_new, snow, shumdiag_perma, &
747	soilcap, soilflx, depth_organic_soil, &
748	stempdiag, gtemp, &
749	mc_layh, mcl_layh, soilmoist, njsc , &
750	frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
751	snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb, &
752	som_total, &
753	veget_max, mc_layh_pft, mcl_layh_pft, soilmoist_pft)
754
755
756	!! 1.12 Initialize river routing
757	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
758	!! 1.12.1 Initialize river routing
759	CALL routing_initialize( kjit, kjpindex, index, &
760	rest_id, hist_id, hist2_id, lalo, &
761	neighbours, resolution, contfrac, stempdiag, &
762	returnflow, reinfiltration, irrigation, riverflow, &
763	coastalflow, flood_frac, flood_res )
764	ELSE
765	!! 1.12.2 No routing, set variables to zero
766	riverflow(:) = zero
767	coastalflow(:) = zero
768	returnflow(:) = zero
769	reinfiltration(:) = zero
770	irrigation(:) = zero
771	flood_frac(:) = zero
772	flood_res(:) = zero
773
774	CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba)
775	CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba)
776
777	ENDIF
778
779	!! 1.13 Write internal variables to output fields
780	z0m_out(:) = z0m(:)
781	z0h_out(:) = z0h(:)
782	emis_out(:) = emis(:)
783	qsurf_out(:) = qsurf(:)
784
785	!! 2. Output variables only once
786	zmaxh_glo(:) = zmaxh
787	CALL xios_orchidee_send_field("zmaxh",zmaxh_glo)
788
789	IF (printlev_loc>=3) WRITE(numout,*) 'sechiba_initialize done'
790
791	END SUBROUTINE sechiba_initialize
792
793	!! ==============================================================================================================================\n
794	!! SUBROUTINE : sechiba_main
795	!!
796	!>\BRIEF Main routine for the sechiba module performing three functions:
797	!! calculating temporal evolution of all variables and preparation of output and
798	!! restart files (during the last call only)
799	!!
800	!!\n DESCRIPTION : Main routine for the sechiba module.
801	!! One time step evolution consists of:
802	!! - call sechiba_var_init to do some initialization,
803	!! - call slowproc_main to do some daily calculations
804	!! - call diffuco_main for diffusion coefficient calculation,
805	!! - call enerbil_main for energy budget calculation,
806	!! - call hydrol_main for hydrologic processes calculation,
807	!! - call condveg_main for surface conditions such as roughness, albedo, and emmisivity,
808	!! - call thermosoil_main for soil thermodynamic calculation,
809	!! - call sechiba_end to swap previous to new fields.
810	!!
811	!! RECENT CHANGE(S): None
812	!!
813	!! MAIN OUTPUT VARIABLE(S): Hydrological variables (:: coastalflow and :: riverflow),
814	!! components of the energy budget (:: tsol_rad, :: vevapp, :: fluxsens,
815	!! :: temp_sol_new and :: fluxlat), surface characteristics (:: z0_out, :: emis_out,
816	!! :: tq_cdrag and :: albedo) and land use related CO2 fluxes (:: netco2flux and
817	!! :: fco2_lu, :: fco2_wh, ::fco2_ha)
818	!!
819	!! REFERENCE(S) :
820	!!
821	!! FLOWCHART :
822	!! \latexonly
823	!! \includegraphics[scale = 0.5]{sechibamainflow.png}
824	!! \endlatexonly
825	!! \n
826	!_ ================================================================================================================================
827
828	SUBROUTINE sechiba_main (kjit, kjpij, kjpindex, index, &
829	& ldrestart_read, ldrestart_write, &
830	& lalo, contfrac, neighbours, resolution,&
831	& zlev, u, v, qair, temp_air, epot_air, ccanopy, &
832	& tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
833	& precip_rain, precip_snow, lwdown, swnet, swdown, coszang, pb, &
834	& vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
835	& netco2flux, fco2_lu, fco2_wh, fco2_ha, &
836	& tsol_rad, temp_sol_new, qsurf_out, albedo, emis_out, z0m_out, z0h_out,&
837	& veget_out, lai_out, height_out, &
838	& rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
839
840	!! 0.1 Input variables
841
842	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
843	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
844	!! (unitless)
845	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
846	!! (unitless)
847	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
848	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (unitless)
849	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier (unitless)
850	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
851	!! (unitless)
852	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _History_ file identifier
853	!! (unitless)
854	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file
855	!! identifier (unitless)
856	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
857	!! (true/false)
858	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
859	!! (true/false)
860	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographic coordinates (latitude,longitude)
861	!! for grid cells (degrees)
862	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid
863	!! (unitless, 0-1)
864	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the pixels on the map.
865	!! Sechiba uses a reduced grid excluding oceans
866	!! ::index contains the indices of the
867	!! terrestrial pixels only! (unitless)
868	INTEGER(i_std), DIMENSION(kjpindex,NbNeighb), INTENT(in) :: neighbours !! Neighboring grid points if land!(unitless)
869	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size in x and y of the grid (m)
870
871	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed in direction u
872	!! @tex $(m.s^{-1})$ @endtex
873	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed in direction v
874	!! @tex $(m.s^{-1})$ @endtex
875	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
876	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
877	!! @tex $(kg kg^{-1})$ @endtex
878	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
879	!! @tex $(kg m^{-2})$ @endtex
880	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
881	!! @tex $(kg m^{-2})$ @endtex
882	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux
883	!! @tex $(W m^{-2})$ @endtex
884	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: coszang !! Cosine of the solar zenith angle (unitless)
885	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux
886	!! @tex $(W m^{-2})$ @endtex
887	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux
888	!! @tex $(W m^{-2})$ @endtex
889	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
890	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (??J)
891	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration in the canopy (ppm)
892	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! Coefficients A for T from the Planetary
893	!! Boundary Layer
894	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! Coefficients A for q from the Planetary
895	!! Boundary Layer
896	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! Coefficients B for T from the Planetary
897	!! Boundary Layer
898	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! Coefficients B for q from the Planetary
899	!! Boundary Layer
900	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
901
902
903	!! 0.2 Output variables
904
905	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: coastalflow !! Outflow on coastal points by small basins.
906	!! This is the water which flows in a disperse
907	!! way into the ocean
908	!! @tex $(kg dt_routing^{-1})$ @endtex
909	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: riverflow !! Outflow of the major rivers.
910	!! The flux will be located on the continental
911	!! grid but this should be a coastal point
912	!! @tex $(kg dt_routing^{-1})$ @endtex
913	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature
914	!! @tex $(W m^{-2})$ @endtex
915	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation
916	!! @tex $(kg m^{-2} days^{-1})$ @endtex
917
918	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf_out !! Surface specific humidity
919	!! @tex $(kg kg^{-1})$ @endtex
920	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m_out !! Surface roughness momentum (output diagnostic, m)
921	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h_out !! Surface roughness heat (output diagnostic, m)
922	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Surface albedo for visible and near-infrared
923	!! (unitless, 0-1)
924	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux
925	!! @tex $(W m^{-2})$ @endtex
926	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux
927	!! @tex $(W m^{-2})$ @endtex
928	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis_out !! Emissivity (unitless)
929	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: netco2flux !! Sum CO2 flux over PFTs
930	!! (gC/m2/dt_sechiba)
931	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_lu !! Land Cover Change CO2 flux (gC/m2/one_day)
932	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_wh !! Wood harvest CO2 flux (gC/m2/one_day)
933	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_ha !! Crop harvest CO2 flux (gC/m2/one_day)
934	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget_out !! Fraction of vegetation type (unitless, 0-1)
935	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: lai_out !! Leaf area index (m^2 m^{-2})
936	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height_out !! Vegetation Height (m)
937	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New ground temperature (K)
938
939	!! 0.3 Modified
940
941	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: tq_cdrag !! Surface drag coefficient (-)
942
943	!! 0.4 local variables
944
945	INTEGER(i_std) :: ji, jv, ilevel !! Index (unitless)
946	REAL(r_std), DIMENSION(kjpindex) :: histvar !! Computations for history files (unitless)
947	REAL(r_std), DIMENSION(kjpindex,nlut) :: histvar2 !! Computations for history files (unitless)
948	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
949	REAL(r_std), DIMENSION(kjpindex) :: sum_treefrac !! Total fraction occupied by trees (0-1, uniless)
950	REAL(r_std), DIMENSION(kjpindex) :: sum_grassfracC3 !! Total fraction occupied by C3 grasses (0-1, unitless)
951	REAL(r_std), DIMENSION(kjpindex) :: sum_grassfracC4 !! Total fraction occupied by C4 grasses (0-1, unitless)
952	REAL(r_std), DIMENSION(kjpindex) :: sum_cropfracC3 !! Total fraction occupied by C3 crops (0-1, unitess)
953	REAL(r_std), DIMENSION(kjpindex) :: sum_cropfracC4 !! Total fraction occupied by C4 crops (0-1, unitess)
954	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracNdlEvg!! Total fraction occupied by treeFracNdlEvg (0-1, unitess)
955	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracBdlEvg!! Total fraction occupied by treeFracBdlEvg (0-1, unitess)
956	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracNdlDcd!! Total fraction occupied by treeFracNdlDcd (0-1, unitess)
957	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracBdlDcd!! Total fraction occupied by treeFracBdlDcd (0-1, unitess)
958	REAL(r_std),DIMENSION (kjpindex,nvm) :: tmc_pft !! Total soil water per PFT (mm/m2)
959	REAL(r_std),DIMENSION (kjpindex,nvm) :: drainage_pft !! Drainage per PFT (mm/m2)
960	REAL(r_std),DIMENSION (kjpindex,nvm) :: runoff_pft !! Runoff per PFT (mm/m2)
961	REAL(r_std),DIMENSION (kjpindex,nvm) :: swc_pft !! Relative Soil water content [tmcr:tmcs] per pft (-)
962	REAL(r_std), DIMENSION(kjpindex) :: grndflux !! Net energy into soil (W/m2)
963	REAL(r_std), DIMENSION(kjpindex,nsnow) :: snowliq !! Liquid water content (m)
964	!+++CHECK+++
965	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
966	! when running a larger domain. Needs to be corrected when implementing
967	! a global use of the multi-layer energy budget.
968	!!$ REAL(r_std),DIMENSION(nlevels_tot) :: Light_Abs_Tot_mean!! total light absorption for a given canopy level
969	!!$ REAL(r_std),DIMENSION(nlevels_tot) :: Light_Alb_Tot_mean!! total albedo for a given level
970	!+++++++++++
971	REAL(r_std), DIMENSION(kjpindex) :: snow_age_diag !! Only for diag, contains xios_default_val
972	REAL(r_std), DIMENSION(kjpindex,nnobio) :: snow_nobio_age_diag !! Only for diag, contains xios_default_val
973	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
974	REAL(r_std), DIMENSION(kjpindex,nlut) :: gpplut !! GPP on landuse tile, only for diagnostics
975	REAL(r_std), DIMENSION(kjpindex,nslm,nvm) :: mc_layh_pft !! As mc_layh per pft
976	REAL(r_std), DIMENSION(kjpindex,nslm,nvm) :: mcl_layh_pft !! As mcl_layh per pft
977	REAL(r_std), DIMENSION(kjpindex,nslm,nvm) :: soilmoist_pft !! As soilmoist per pft
978	REAL(r_std), DIMENSION(kjpindex) :: vbeta2sum !! sum of vbeta2 coefficients across all PFTs (-)
979	REAL(r_std), DIMENSION(kjpindex) :: vbeta3sum !! sum of vbeta3 coefficients across all PFTs (-)
980	REAL(r_std), DIMENSION(kjpindex) :: netrad !! Net radiation (W m^{-2})
981	REAL(r_std), DIMENSION(kjpindex) :: lwabs !! LW radiation absorbed by the surface (W m^{-2})
982	REAL(r_std), DIMENSION(kjpindex) :: lwnet !! Net Long-wave radiation (W m^{-2})
983	REAL(r_std), DIMENSION(kjpindex) :: fluxsubli !! Energy of sublimation (mm day^{-1})
984
985	REAL(r_std),DIMENSION (kjpindex,nvm) :: vbeta23 !! Beta for fraction of wetted foliage that will
986	!! transpire once intercepted water has evaporated (-)
987	REAL(r_std), DIMENSION (kjpindex,nvm,nlai) :: leaf_ci !! intercellular CO2 concentration (ppm)
988	REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot) :: gs_distribution
989	REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot) :: gs_diffuco_output
990	REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot) :: gstot_component
991	REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot) :: gstot_frac
992	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_vbeta3
993	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_rveget
994	INTEGER :: ilev,ivm,ipts
995	LOGICAL, DIMENSION(kjpindex,nvm) :: failed_vegfrac !! Pixels and pfts for which the model failed to
996	!! find a PFT were some residual fraction could be added (true/false)
997	REAL(r_std),DIMENSION (kjpindex,nvm) :: cresist !!coefficient for resistances (??)
998	REAL(r_std), DIMENSION(kjpindex) :: mcs_hydrol !! Saturated volumetric water content output to be used in stomate_soilcarbon
999	REAL(r_std), DIMENSION(kjpindex) :: mcfc_hydrol !! Volumetric water content at field capacity output to be used in stomate_soilcarbon
1000
1001	REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot+1) :: info_limitphoto !!Save information of limitation for photosynthesis
1002	REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot) :: JJ_out
1003	REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot) :: assimi_lev
1004	LOGICAL :: hydrol_flag !! flag that 'trips' the energy budget for each grid square,
1005	!! but should be false when mleb is called from sechiba
1006	LOGICAL, DIMENSION (kjpindex) :: hydrol_flag2 !! flag that 'trips' the energy budget for each grid square,
1007	!! hydrol_flag2 on the full global domain for all processors
1008	!! but should be false when mleb is called from sechiba
1009	LOGICAL, DIMENSION (kjpindex,nvm) :: hydrol_flag3 !! flag that 'trips' the energy budget for each grid square and PFT
1010	!! but should be false when mleb is called from sechiba
1011
1012	REAL(r_std), DIMENSION (kjpindex,nvm) :: temporary_array !! A temporary array
1013	REAL(r_std), DIMENSION (kjpindex) :: zenith_angle
1014	REAL(r_std),DIMENSION (kjpindex,nvm) :: co2_flux !! CO2 flux per average ground area (gC/m2/one_day)
1015
1016	!_ ================================================================================================================================
1017
1018	IF (printlev_loc>=3) WRITE(numout,*) 'Start sechiba_main kjpindex =',kjpindex
1019
1020	IF (FirstTsDay) THEN
1021	CALL solarang_noon(kjpindex, lalo, coszang_noon)
1022	END IF
1023	zenith_angle(:) = ACOS(coszang_noon)*(180/pi)
1024	WHERE (zenith_angle(:).LE.zero)
1025	zenith_angle(:) = xios_default_val
1026	ENDWHERE
1027	CALL xios_orchidee_send_field("zenith_ang_noon",zenith_angle(:))
1028
1029	!! 1. Initialize variables at each time step
1030	CALL sechiba_var_init (kjpindex, rau, pb, temp_air)
1031
1032	!! 2. Compute diffusion coefficients
1033	CALL diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, &
1034	& zlev, z0m, z0h, roughheight, temp_sol, temp_air, temp_growth, rau, tq_cdrag, qsurf, qair, pb , &
1035	& evap_bare_lim, evap_bare_lim_ns, evapot, evapot_corr, snow, flood_frac, flood_res, &
1036	& frac_nobio, snow_nobio, totfrac_nobio, &
1037	& swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai, &
1038	& qsintveg, qsintmax, assim_param, &
1039	& vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, gsmean, rveget, rstruct, cimean, gpp, cresist, &
1040	& lalo, neighbours, resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
1041	& hist_id, hist2_id, vbeta2sum, vbeta3sum, Light_Abs_Tot, &
1042	& Light_Tran_Tot, lai_per_level, vbeta23, leaf_ci, &
1043	& gs_distribution, gs_diffuco_output, gstot_component, gstot_frac,&
1044	& warnings, u_speed, profile_vbeta3, profile_rveget, &
1045	& delta_c13_assim, leaf_ci_out, info_limitphoto, JJ_out, assimi_lev)
1046
1047	IF ( ok_c13 ) THEN
1048
1049	! When there is no photosynthesis, carbon isotopic values is not calulated.
1050	! Because daily mean of carbon isotopic value in ORCHIDEE contained night values,
1051	! to avoide biased daily mean value by zeros, daily mean of carbon isotopic
1052	! discriminations should be calculated only when there is photosynthesis
1053
1054	gpp_day(:,:) = zero
1055
1056	DO ipts = 1,kjpindex
1057	DO jv = 1,nvm
1058	IF (gpp(ipts,jv) .LT. min_sechiba) THEN
1059	gpp_day(ipts,jv) = zero
1060	ELSE
1061	gpp_day(ipts,jv) = gpp_day(ipts,jv) + 1
1062	END IF
1063	END DO
1064	END DO
1065
1066	END IF ! ok_c13
1067
1068	!! 3. Compute energy balance
1069	IF (ok_hydrol_arch) THEN
1070
1071	CALL hydraulic_arch_main( &
1072	kjit, kjpij, kjpindex, rest_id, &
1073	rest_id_stom, hist_id, hist2_id, hist_id_stom, &
1074	hist_id_stom_IPCC, index, indexveg, ldrestart_read, &
1075	ldrestart_write, njsc, u, v, &
1076	qair, precip_rain, precip_snow, lwdown, &
1077	swnet, swdown, temp_air, &
1078	petAcoef, peqAcoef, petBcoef, peqBcoef, &
1079	pb, tq_cdrag, epot_air, coszang, &
1080	ccanopy, emis, soilcap, soilflx, &
1081	rau, temp_growth, vbeta23, humrel, &
1082	lai, max_height_store, gs_distribution, gs_diffuco_output, &
1083	gstot_component, gstot_frac, assim_param, circ_class_biomass, &
1084	circ_class_n, ksoil, lai_per_level, laieff_isotrop, &
1085	Light_Abs_Tot, Light_Tran_Tot, qsintmax, qsintveg, &
1086	snowdz, stempdiag, swc, &
1087	vbeta1, vbeta2, vbeta4, vbeta5, &
1088	veget_max, veget, z_array_out, &
1089	temp_sol_new, tsol_rad, temp_sol_add, vevapflo, &
1090	vevapnu, vevapsno, stressed, unstressed, &
1091	transpot, vevapwet, transpir, &
1092	e_frac, &
1093	vbeta2sum, vbeta3sum, fluxsens, fluxlat, &
1094	vevapp, vbeta, evapot, evapot_corr, &
1095	qsurf, temp_sol, pgflux, u_speed, &
1096	vbeta3, vbeta3pot, rveget, rstruct, &
1097	cimean, gsmean, gpp, transpir_supply, &
1098	leaf_ci, warnings, JJ_out, assimi_lev, &
1099	cresist, info_limitphoto, vessel_loss, root_profile)
1100
1101	ELSEIF(.NOT. ok_hydrol_arch .AND. .NOT. ok_mleb) THEN
1102
1103	CALL enerbil_main (kjit, kjpindex, netrad, lwabs, lwnet, fluxsubli, index, &
1104	indexveg, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, qair, &
1105	peqAcoef, peqBcoef, pb, rau, vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, &
1106	vbeta5, emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
1107	vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo, temp_sol, tsol_rad, &
1108	temp_sol_new, qsurf, evapot, evapot_corr, precip_rain, pgflux, &
1109	snowdz, temp_sol_add, qair_new, qsol_sat_new)
1110
1111	! Write enerbil output. Because enerbil_main is called twice in the case
1112	! with water stress when ok_hydrol_arch, calling this module in the
1113	! enrebil_main can cause a problem. For now, the module is placed separately.
1114	CALL enerbil_write (kjit, kjpindex, hist_id, hist2_id, index, &
1115	evapot, evapot_corr, fluxsubli, lwabs, lwnet, netrad, transpir, &
1116	vevapflo, vevapnu, vevapp, vevapsno, vevapwet)
1117
1118	ELSEIF(.NOT. ok_hydrol_arch .AND. ok_mleb) THEN
1119
1120	! These hydrol flag will always be false, when the model is not run with
1121	! the hydraulic architecture.
1122	hydrol_flag = .FALSE.
1123	hydrol_flag2 = .FALSE.
1124	hydrol_flag3 = .FALSE.
1125	!+++CHECK+++
1126	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
1127	! when running a larger domain. Needs to be corrected when implementing
1128	! a global use of the multi-layer energy budget.
1129	CALL mleb_main (kjit, kjpindex, ldrestart_read, ldrestart_write, &
1130	& index, indexveg, lwdown, swnet, swdown, epot_air, temp_air, &
1131	& u, v, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, pb, rau, &
1132	& vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
1133	& emis, soilflx, soilcap, tq_cdrag, humrel, &
1134	& fluxsens, fluxlat, vevapp, transpir, transpot, vevapnu,vevapwet, &
1135	& vevapsno, vevapflo, temp_sol, tsol_rad, temp_sol_new, qsurf, &
1136	& evapot, evapot_corr, rest_id, hist_id, hist2_id, &
1137	& ok_mleb_history_file, &
1138	& transpir_supply, hydrol_flag, hydrol_flag2, hydrol_flag3,vbeta2sum, &
1139	& vbeta3sum, veget_max, qsol_sat_new, qair_new, &
1140	& veget, &
1141	!!$ ! Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
1142	laieff_isotrop,&
1143	& z_array_out, transpir_supply_column, u_speed, &
1144	& profile_vbeta3, profile_rveget, max_height_store, &
1145	& precip_rain, pgflux, snowdz, temp_sol_add)
1146	!+++++++++++
1147
1148
1149	ELSE
1150
1151	CALL ipslerr_p(3,'sechiba_main','This should not happen',&
1152	'Please chech the configurations for the eneryg',&
1153	'We need energy budget calculations')
1154
1155
1156	END IF
1157
1158	!! 4. Compute hydrology
1159	!! 4.1 Water balance from CWRR module (11 soil layers)
1160	CALL hydrol_main (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjit, kjpindex, &
1161	& index, indexveg, indexsoil, indexlayer, indexnslm, &
1162	& temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, &
1163	& qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age, &
1164	& tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, &
1165	& humrel, vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, evap_bare_lim_ns, flood_frac, flood_res, &
1166	& shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope,&
1167	& rest_id, hist_id, hist2_id,&
1168	& contfrac, stempdiag, &
1169	& temp_air, pb, u, v, tq_cdrag, swnet, pgflux, &
1170	& snowrho,snowtemp,snowgrain,snowdz,snowheat,snowliq, &
1171	& grndflux,gtemp,tot_bare_soil, &
1172	& lambda_snow,cgrnd_snow,dgrnd_snow,frac_snow_veg,temp_sol_add, &
1173	& mc_layh, mcl_layh, tmc_pft, drainage_pft, runoff_pft, swc_pft, soilmoist, &
1174	& mc_layh_s, mcl_layh_s, soilmoist_s, swc, e_frac, ksoil, &
1175	& mcs_hydrol, mcfc_hydrol, altmax, root_profile, root_depth, &
1176	& circ_class_biomass, us)
1177
1178
1179	!! 6. Compute surface variables (emissivity, albedo and roughness)
1180	!+++CHECK+++
1181	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
1182	! when running a larger domain. Needs to be corrected when implementing
1183	! a global use of the multi-layer energy budget.
1184	CALL condveg_main (kjit, kjpindex, index, rest_id, hist_id, hist2_id, &
1185	lalo, neighbours, resolution, contfrac, &
1186	veget, veget_max, frac_nobio, totfrac_nobio, &
1187	zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
1188	drysoil_frac, height, height_dom, snowdz, snowrho, tot_bare_soil, &
1189	temp_air, pb, u, v, &
1190	lai, &
1191	emis, albedo, z0m, z0h, roughheight, &
1192	frac_snow_veg, frac_snow_nobio, coszang, &
1193	Light_Abs_Tot, Light_Tran_Tot, laieff_fit, &
1194	!!$ ! Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
1195	laieff_isotrop)
1196	!+++++++++++
1197
1198	!! 7. Compute soil thermodynamics
1199	DO jv = 1,nvm
1200	mc_layh_pft(:,:,jv) = mc_layh_s(:,:,pref_soil_veg(jv))
1201	mcl_layh_pft(:,:,jv) = mcl_layh_s(:,:,pref_soil_veg(jv))
1202	soilmoist_pft(:,:,jv) = soilmoist_s(:,:,pref_soil_veg(jv))
1203	END DO
1204
1205	CALL thermosoil_main (kjit, kjpindex, &
1206	index, indexgrnd, &
1207	temp_sol_new, snow, soilcap, soilflx, &
1208	shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
1209	snowdz,snowrho,snowtemp,gtemp,pb,&
1210	mc_layh, mcl_layh, soilmoist, &
1211	mc_layh_pft, mcl_layh_pft,soilmoist_pft, njsc, &
1212	depth_organic_soil, heat_Zimov, tdeep, hsdeep,&
1213	som_total, veget_max, &
1214	frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
1215	lambda_snow, cgrnd_snow, dgrnd_snow)
1216
1217
1218	!! 8. Compute river routing
1219	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
1220	!! 8.1 River routing
1221	CALL routing_main (kjit, kjpindex, index, &
1222	& lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
1223	& drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, &
1224	& stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
1225	ELSE
1226	!! 8.2 No routing, set variables to zero
1227	riverflow(:) = zero
1228	coastalflow(:) = zero
1229	returnflow(:) = zero
1230	reinfiltration(:) = zero
1231	irrigation(:) = zero
1232	flood_frac(:) = zero
1233	flood_res(:) = zero
1234
1235	CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba)
1236	CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba)
1237	ENDIF
1238
1239
1240	!! 9. Compute slow processes (i.e. 'daily' and annual time step)
1241	CALL slowproc_main (kjit, kjpij, kjpindex, njsc, &
1242	index, indexveg, lalo, neighbours, resolution, contfrac, soiltile, fraclut, nwdFraclut, &
1243	temp_air, temp_sol, stempdiag, &
1244	vegstress, humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, pb, gpp, &
1245	tmc_pft, drainage_pft, runoff_pft, swc_pft, deadleaf_cover, &
1246	assim_param, qsintveg, &
1247	frac_age, height, lai, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
1248	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
1249	co2_flux, fco2_lu, fco2_wh, fco2_ha, temp_growth, tot_bare_soil, &
1250	tdeep, hsdeep, snow, heat_Zimov, &
1251	sfluxCH4_deep, sfluxCO2_deep, &
1252	som_total, snowdz, snowrho, altmax, depth_organic_soil, &
1253	circ_class_biomass, circ_class_n, &
1254	lai_per_level, max_height_store, laieff_fit, laieff_isotrop, &
1255	z_array_out, transpir, transpir_mod, &
1256	transpir_supply, vir_transpir_supply, coszang, coszang_noon, &
1257	stressed, unstressed, Light_Tran_Tot, &
1258	u, v, loss_gain, veget_max_new, frac_nobio_new, failed_vegfrac, &
1259	mcs_hydrol, mcfc_hydrol,vessel_loss, height_dom, root_profile, &
1260	root_depth, us, Pgap_cumul)
1261
1262	!+++CHECK+++
1263	! slowproc_main basically ends with a CALL to slowproc_veget and
1264	! check_veget if do_slow. LCC takes place at the end of the first day
1265	! so one could expect do_slow = TRUE. Are the following CALLs really
1266	! needed. Seems that they should not do too much. Consider deleting.
1267
1268	! Update variables such as veget, frac_nobio, soiltile, lai_per_level
1269	! after LCC has been done in stomatelpj. Some of these variables are
1270	! used in calculating variables for the history files.
1271	IF (done_stomate_lcchange) THEN
1272
1273	!! Calculation of veget and soiltile.
1274	CALL slowproc_veget (kjpindex, lai_per_level, z_array_out, &
1275	coszang_noon, circ_class_biomass, circ_class_n, frac_nobio, totfrac_nobio, &
1276	veget_max, veget, soiltile, tot_bare_soil, &
1277	fraclut, nwdFraclut, Pgap_cumul)
1278
1279	!! Do some basic tests on the surface fractions updated above
1280	CALL check_veget(kjpindex, frac_nobio, veget_max, &
1281	veget, tot_bare_soil, soiltile, failed_vegfrac)
1282
1283	done_stomate_lcchange = .FALSE.
1284
1285	END IF
1286	!+++++++++++
1287
1288	!! Compute global CO2 flux
1289	! This value has to be calculated every half hour but co2_flux is calculated
1290	! in stomate_lpj and therefore only defined once per day. To avoid
1291	! that it is zero for the first 47 time steps of a day and gets a value
1292	! only the last (48th) time step this value is send to the restart files.
1293	! Expressed in gC m-2 day-1 and send to orchideedriver.f90 in this unit.
1294	! Because of this approach netco2flux and c02_flux are 23.5 hours out
1295	! of sink. co2_flux states that the co2_flux for that day is given by
1296	! the variable co2_flux. For the same day the variable netco2flux uses the
1297	! 47 times the co2_flux value from the day before. Only at time step
1298	! 48 the value of the actual day is used. Not is not a problem but it
1299	! is good to realize when plotting and comparing the values of co2_flux
1300	! (given as nee in the sechiba history files) and netco2flux. Given as
1301	! (netco2flux in the sechiba history files).
1302	!+++CHECK+++
1303	! co2_flux is also used to calculate the mass balance closure in
1304	! stomate_lpj and therefore contains the atm_to_bm fluxes. These fluxes
1305	! are substantial for croplands. The co2_flux in ORC2.2 does not
1306	! atm_to_bm. For that reason the C-cycle im LMDzOR using that version
1307	! of ORCHIDEE cannot be closed. With this approach the C-cycle is closed
1308	! on the ORCHIDEE side but maybe LMDz does not like the sudden spikes
1309	! caused by atm_to_bm.
1310	! co2_flux does not contain any information on harvest and products
1311	! even if this is the variable that is needed in LMDzOR it should not be
1312	! called net co2 flux because it is not at the all the net flux. I [SL]
1313	! strongly doubts how the current definition of netco2flux could be
1314	! useful in LMDzOR.
1315	netco2flux(:) = zero
1316	DO jv = 2,nvm
1317	netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*(1-totfrac_nobio)
1318	ENDDO
1319	!+++++++++++
1320
1321	! Write to sechiba history file for debugging purpose
1322	! convert the units from gC m-2 day-1 to kgC m-2 day-1 as for the other
1323	! C-fluxes in sechiba_history.
1324	CALL xios_orchidee_send_field("netco2flux",netco2flux(:)/1.e3/one_day)
1325
1326
1327	!! 10. Update the temperature (temp_sol) with newly computed values
1328	CALL sechiba_end (kjpindex, temp_sol_new, temp_sol)
1329
1330	!! 14. If it is the last time step, write restart files
1331	IF (ldrestart_write) THEN
1332	!+++CHECK+++
1333	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
1334	! when running a larger domain. Needs to be corrected when implementing
1335	! a global use of the multi-layer energy budget.
1336	CALL sechiba_finalize( &
1337	kjit, kjpij, kjpindex, index, rest_id, &
1338	tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad, &
1339	albedo &
1340	!!$ , Light_Abs_Tot_mean, Light_Alb_Tot_mean &
1341	)
1342	!+++++++++++
1343	END IF
1344
1345	!! 11. Write internal variables to output fields
1346	z0m_out(:) = z0m(:)
1347	z0h_out(:) = z0h(:)
1348	emis_out(:) = emis(:)
1349	qsurf_out(:) = qsurf(:)
1350	veget_out(:,:) = veget(:,:)
1351	lai_out(:,:) = lai(:,:)
1352	height_out(:,:) = height(:,:)
1353
1354	!! 12. Write global variables to history files
1355	sum_treefrac(:) = zero
1356	sum_grassfracC3(:) = zero
1357	sum_grassfracC4(:) = zero
1358	sum_cropfracC3(:) = zero
1359	sum_cropfracC4(:) = zero
1360	sum_treeFracNdlEvg(:) = zero
1361	sum_treeFracBdlEvg(:) = zero
1362	sum_treeFracNdlDcd(:) = zero
1363	sum_treeFracBdlDcd(:) = zero
1364	DO jv = 2, nvm
1365	IF (is_tree(jv) .AND. natural(jv)) THEN
1366	sum_treefrac(:) = sum_treefrac(:) + veget_max(:,jv)
1367	ELSE IF ((.NOT. is_tree(jv)) .AND. natural(jv)) THEN
1368	! Grass
1369	IF (is_c4(jv)) THEN
1370	sum_grassfracC4(:) = sum_grassfracC4(:) + veget_max(:,jv)
1371	ELSE
1372	sum_grassfracC3(:) = sum_grassfracC3(:) + veget_max(:,jv)
1373	END IF
1374	ELSE
1375	! Crop and trees not natural
1376	IF (is_c4(jv)) THEN
1377	sum_cropfracC4(:) = sum_cropfracC4(:) + veget_max(:,jv)
1378	ELSE
1379	sum_cropfracC3(:) = sum_cropfracC3(:) + veget_max(:,jv)
1380	END IF
1381	ENDIF
1382
1383	IF (is_tree(jv)) THEN
1384	IF (is_evergreen(jv)) THEN
1385	IF (is_needleleaf(jv)) THEN
1386	! Fraction for needleleaf evergreen trees (treeFracNdlEvg)
1387	sum_treeFracNdlEvg(:) = sum_treeFracNdlEvg(:) + veget_max(:,jv)
1388	ELSE
1389	! Fraction for broadleaf evergreen trees (treeFracBdlEvg)
1390	sum_treeFracBdlEvg(:) = sum_treeFracBdlEvg(:) + veget_max(:,jv)
1391	END IF
1392	ELSE IF (is_deciduous(jv)) THEN
1393	IF (is_needleleaf(jv)) THEN
1394	! Fraction for needleleaf deciduous trees (treeFracNdlDcd)
1395	sum_treeFracNdlDcd(:) = sum_treeFracNdlDcd(:) + veget_max(:,jv)
1396	ELSE
1397	! Fraction for broadleafs deciduous trees (treeFracBdlDcd)
1398	sum_treeFracBdlDcd(:) = sum_treeFracBdlDcd(:) + veget_max(:,jv)
1399	END IF
1400	END IF
1401	END IF
1402	ENDDO
1403
1404	histvar(:)=zero
1405	DO jv = 2, nvm
1406	IF (is_deciduous(jv)) THEN
1407	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
1408	ENDIF
1409	ENDDO
1410	CALL xios_orchidee_send_field("treeFracPrimDec",histvar)
1411
1412	histvar(:)=zero
1413	DO jv = 2, nvm
1414	IF (is_evergreen(jv)) THEN
1415	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
1416	ENDIF
1417	ENDDO
1418	CALL xios_orchidee_send_field("treeFracPrimEver",histvar)
1419
1420	histvar(:)=zero
1421	DO jv = 2, nvm
1422	IF ( .NOT.(is_c4(jv)) ) THEN
1423	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
1424	ENDIF
1425	ENDDO
1426	CALL xios_orchidee_send_field("c3PftFrac",histvar)
1427
1428	histvar(:)=zero
1429	DO jv = 2, nvm
1430	IF ( is_c4(jv) ) THEN
1431	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
1432	ENDIF
1433	ENDDO
1434	CALL xios_orchidee_send_field("c4PftFrac",histvar)
1435
1436	CALL xios_orchidee_send_field("temp_sol_new",temp_sol_new)
1437	CALL xios_orchidee_send_field("fluxsens",fluxsens)
1438	CALL xios_orchidee_send_field("fluxlat",fluxlat)
1439
1440
1441	! Add XIOS default value where no snow
1442	DO ji=1,kjpindex
1443	IF (snow(ji) .GT. zero) THEN
1444	snow_age_diag(ji) = snow_age(ji)
1445	snow_nobio_age_diag(ji,:) = snow_nobio_age(ji,:)
1446
1447	snowage_glob(ji) = snow_age(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji)) + &
1448	SUM(snow_nobio_age(ji,:)frac_snow_nobio(ji,:)frac_nobio(ji,:))
1449	IF (snowage_glob(ji) .NE. 0) snowage_glob(ji) = snowage_glob(ji) / &
1450	(frac_snow_veg(ji)(1-totfrac_nobio(ji)) + SUM(frac_snow_nobio(ji,:)frac_nobio(ji,:)))
1451	ELSE
1452	snow_age_diag(ji) = xios_default_val
1453	snow_nobio_age_diag(ji,:) = xios_default_val
1454	snowage_glob(ji) = xios_default_val
1455	END IF
1456	END DO
1457
1458	CALL xios_orchidee_send_field("snow",snow)
1459	CALL xios_orchidee_send_field("snowage",snow_age_diag)
1460	CALL xios_orchidee_send_field("snownobio",snow_nobio)
1461	CALL xios_orchidee_send_field("snownobioage",snow_nobio_age_diag)
1462	CALL xios_orchidee_send_field("snowage_glob",snowage_glob)
1463
1464	CALL xios_orchidee_send_field("frac_snow", SUM(frac_snow_nobio,2)totfrac_nobio+frac_snow_veg(1-totfrac_nobio))
1465	CALL xios_orchidee_send_field("frac_snow_veg", frac_snow_veg)
1466	CALL xios_orchidee_send_field("frac_snow_nobio", frac_snow_nobio)
1467	CALL xios_orchidee_send_field("frac_snow", SUM(frac_snow_nobio,2)totfrac_nobio+frac_snow_veg(1-totfrac_nobio))
1468	CALL xios_orchidee_send_field("frac_snow_veg", frac_snow_veg)
1469	CALL xios_orchidee_send_field("frac_snow_nobio", frac_snow_nobio)
1470	CALL xios_orchidee_send_field("pgflux",pgflux)
1471	CALL xios_orchidee_send_field("reinf_slope",reinf_slope)
1472	CALL xios_orchidee_send_field("njsc",REAL(njsc, r_std))
1473	CALL xios_orchidee_send_field("vegetfrac",veget)
1474	CALL xios_orchidee_send_field("maxvegetfrac",veget_max)
1475	CALL xios_orchidee_send_field("nobiofrac",frac_nobio)
1476	CALL xios_orchidee_send_field("soiltile",soiltile)
1477	CALL xios_orchidee_send_field("rstruct",rstruct)
1478
1479	!-
1480	! Write effective LAI to the output file. We choose the isotropic LAI
1481	! here since that is integrated over all solar angles. Comparisons
1482	! come from satellies and are dependent on the viewing angle, which is not
1483	! available from observations, so this becomes the LAI as viewed from
1484	! 60 degrees instead of the typical LAI, which is viewed from 90 degrees.
1485	! laieff_isotrop is calculated for every level of LAI, but we just
1486	! want to print out a single value for the canopy.
1487	! This is a variable that is part stomate, part sechiba. Print it
1488	! out in sechiba because it's used here, even though it depends
1489	! completely on stomate.
1490	temporary_array(:,:)=0.0d0
1491	DO ilevel = 1,nlevels_tot
1492	temporary_array(:,:)=temporary_array(:,:)+laieff_isotrop(:,ilevel,:)
1493	ENDDO
1494	CALL xios_orchidee_send_field("EFFECTIVE_LAI",temporary_array(:,:))
1495
1496	!+++CHECK+++
1497	! AHAAA These values do not account for atm_to_bm. They are incomplete
1498	! and thus wrong!
1499	CALL xios_orchidee_send_field("gpp",gpp/dt_sechiba)
1500	CALL xios_orchidee_send_field("gpp_ipcc2",SUM(gpp,dim=2)/dt_sechiba)
1501
1502	histvar(:)=zero
1503	DO jv = 2, nvm
1504	IF ( .NOT. is_tree(jv) .AND. natural(jv) ) THEN
1505	histvar(:) = histvar(:) + gpp(:,jv)
1506	ENDIF
1507	ENDDO
1508	CALL xios_orchidee_send_field("gppgrass",histvar/dt_sechiba)
1509
1510	histvar(:)=zero
1511	DO jv = 2, nvm
1512	IF ( (.NOT. is_tree(jv)) .AND. (.NOT. natural(jv)) ) THEN
1513	histvar(:) = histvar(:) + gpp(:,jv)
1514	ENDIF
1515	ENDDO
1516	CALL xios_orchidee_send_field("gppcrop",histvar/dt_sechiba)
1517
1518	histvar(:)=zero
1519	DO jv = 2, nvm
1520	IF ( is_tree(jv) ) THEN
1521	histvar(:) = histvar(:) + gpp(:,jv)
1522	ENDIF
1523	ENDDO
1524	CALL xios_orchidee_send_field("gpptree",histvar/dt_sechiba)
1525	!++++++++++++
1526
1527	CALL xios_orchidee_send_field("drysoil_frac",drysoil_frac)
1528	CALL xios_orchidee_send_field("vevapflo",vevapflo/dt_sechiba)
1529	CALL xios_orchidee_send_field("k_litt",k_litt)
1530	CALL xios_orchidee_send_field("beta",vbeta)
1531	CALL xios_orchidee_send_field("vbeta1",vbeta1)
1532	CALL xios_orchidee_send_field("vbeta2",vbeta2)
1533	CALL xios_orchidee_send_field("vbeta3",vbeta3)
1534	CALL xios_orchidee_send_field("vbeta4",vbeta4)
1535	CALL xios_orchidee_send_field("vbeta5",vbeta5)
1536
1537	CALL xios_orchidee_send_field("gsmean",gsmean)
1538	CALL xios_orchidee_send_field("cimean",cimean)
1539	CALL xios_orchidee_send_field("rveget",rveget)
1540
1541	IF ( ok_c13 ) THEN
1542	WHERE (ABS(delta_c13_assim) < min_sechiba)
1543	delta_c13_assim = xios_default_val
1544	ENDWHERE
1545	WHERE (leaf_ci_out < min_sechiba)
1546	leaf_ci_out = xios_default_val
1547	ENDWHERE
1548	WHERE (gpp_day < min_sechiba)
1549	gpp_day = xios_default_val
1550	ENDWHERE
1551	ELSE
1552	delta_c13_assim(:,:)=-9999
1553	leaf_ci_out(:,:)=-9999
1554	gpp_day(:,:)=-9999
1555	ENDIF
1556
1557	CALL xios_orchidee_send_field('delta_c13_assim', delta_c13_assim)
1558	CALL xios_orchidee_send_field('leaf_ci_out', leaf_ci_out)
1559	CALL xios_orchidee_send_field('gpp_day', gpp_day)
1560
1561	histvar(:)=SUM(vevapwet(:,:),dim=2)
1562	CALL xios_orchidee_send_field("evspsblveg",histvar/dt_sechiba)
1563	histvar(:)= vevapnu(:)+vevapsno(:)
1564	CALL xios_orchidee_send_field("evspsblsoi",histvar/dt_sechiba)
1565	histvar(:)=SUM(transpir(:,:),dim=2)
1566	CALL xios_orchidee_send_field("tran",histvar/dt_sechiba)
1567
1568	! For CMIP6 data request: the following fractions are fractions of the total grid-cell,
1569	! which explains the multiplication by contfrac
1570	CALL xios_orchidee_send_field("treeFrac",sum_treefrac(:)100contfrac(:))
1571	CALL xios_orchidee_send_field("grassFracC3",sum_grassfracC3(:)100contfrac(:))
1572	CALL xios_orchidee_send_field("grassFracC4",sum_grassfracC4(:)100contfrac(:))
1573	CALL xios_orchidee_send_field("cropFracC3",sum_cropfracC3(:)100contfrac(:))
1574	CALL xios_orchidee_send_field("cropFracC4",sum_cropfracC4(:)100contfrac(:))
1575	CALL xios_orchidee_send_field("treeFracNdlEvg",sum_treeFracNdlEvg(:)100contfrac(:))
1576	CALL xios_orchidee_send_field("treeFracBdlEvg",sum_treeFracBdlEvg(:)100contfrac(:))
1577	CALL xios_orchidee_send_field("treeFracNdlDcd",sum_treeFracNdlDcd(:)100contfrac(:))
1578	CALL xios_orchidee_send_field("treeFracBdlDcd",sum_treeFracBdlDcd(:)100contfrac(:))
1579
1580	histvar(:)=veget_max(:,1)100contfrac(:)
1581	CALL xios_orchidee_send_field("baresoilFrac",histvar)
1582	histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)100contfrac(:)
1583	CALL xios_orchidee_send_field("residualFrac",histvar)
1584
1585	! For CMIP6 data request: cnc = canopy cover fraction over land area
1586	histvar(:)=zero
1587	DO jv=2,nvm
1588	histvar(:) = histvar(:) + veget_max(:,jv)*100
1589	END DO
1590	CALL xios_orchidee_send_field("cnc",histvar)
1591
1592	CALL xios_orchidee_send_field("tsol_rad",tsol_rad-273.15)
1593	CALL xios_orchidee_send_field("qsurf",qsurf)
1594	CALL xios_orchidee_send_field("emis",emis)
1595	CALL xios_orchidee_send_field("z0m",z0m)
1596	CALL xios_orchidee_send_field("z0h",z0h)
1597	CALL xios_orchidee_send_field("roughheight",roughheight)
1598
1599
1600	! calculate laimean (pixel including no_bio)
1601	histvar(:)=zero
1602	DO ji = 1, kjpindex
1603	IF (SUM(veget_max(ji,:)) > zero) THEN
1604	DO jv=2,nvm
1605	histvar(ji) = histvar(ji) + veget_max(ji,jv)*lai(ji,jv)/ &
1606	SUM(veget_max(ji,:))
1607	END DO
1608	END IF
1609	END DO
1610	CALL xios_orchidee_send_field("lai",lai)
1611	CALL xios_orchidee_send_field("LAImean",histvar)
1612
1613	CALL xios_orchidee_send_field("vevapsno",vevapsno/dt_sechiba)
1614	CALL xios_orchidee_send_field("vevapp",vevapp/dt_sechiba)
1615	CALL xios_orchidee_send_field("vevapnu",vevapnu/dt_sechiba)
1616	CALL xios_orchidee_send_field("transpir",transpir*one_day/dt_sechiba)
1617	CALL xios_orchidee_send_field("inter",vevapwet*one_day/dt_sechiba)
1618
1619	histvar(:)=zero
1620	DO jv=1,nvm
1621	histvar(:) = histvar(:) + vevapwet(:,jv)
1622	ENDDO
1623	CALL xios_orchidee_send_field("ECanop",histvar/dt_sechiba)
1624	histvar(:)=zero
1625	DO jv=1,nvm
1626	histvar(:) = histvar(:) + transpir(:,jv)
1627	ENDDO
1628	CALL xios_orchidee_send_field("TVeg",histvar/dt_sechiba)
1629
1630	!! Calculate diagnostic variables on Landuse tiles for LUMIP/CMIP6
1631
1632	! Calculate fraction of landuse tiles related to the whole grid cell
1633	DO jv=1,nlut
1634	histvar2(:,jv) = fraclut(:,jv) * contfrac(:)
1635	END DO
1636	CALL xios_orchidee_send_field("fraclut",histvar2)
1637
1638	CALL xios_orchidee_send_field("nwdFraclut",nwdFraclut(:,:))
1639
1640	! Calculate GPP on landuse tiles
1641	! val_exp is used as missing value where the values are not known i.e. where the tile is not represented
1642	! or for pasture (id_pst) or urban land (id_urb).
1643	gpplut(:,:)=0
1644	DO jv=1,nvm
1645	IF (natural(jv)) THEN
1646	gpplut(:,id_psl) = gpplut(:,id_psl) + gpp(:,jv)
1647	ELSE
1648	gpplut(:,id_crp) = gpplut(:,id_crp) + gpp(:,jv)
1649	ENDIF
1650	END DO
1651
1652	! Transform from gC/m2/s into kgC/m2/s
1653	WHERE (fraclut(:,id_psl)>min_sechiba)
1654	gpplut(:,id_psl) = gpplut(:,id_psl)/fraclut(:,id_psl)/1000
1655	ELSEWHERE
1656	gpplut(:,id_psl) = xios_default_val
1657	END WHERE
1658	WHERE (fraclut(:,id_crp)>min_sechiba)
1659	gpplut(:,id_crp) = gpplut(:,id_crp)/fraclut(:,id_crp)/1000
1660	ELSEWHERE
1661	gpplut(:,id_crp) = xios_default_val
1662	END WHERE
1663	gpplut(:,id_pst) = xios_default_val
1664	gpplut(:,id_urb) = xios_default_val
1665
1666	CALL xios_orchidee_send_field("gpplut",gpplut)
1667
1668
1669	IF ( .NOT. almaoutput ) THEN
1670	! Write history file in IPSL-format
1671	CALL histwrite_p(hist_id, 'beta', kjit, vbeta, kjpindex, index)
1672	CALL histwrite_p(hist_id, 'z0m', kjit, z0m, kjpindex, index)
1673	CALL histwrite_p(hist_id, 'z0h', kjit, z0h, kjpindex, index)
1674	CALL histwrite_p(hist_id, 'roughheight', kjit, roughheight, kjpindex, index)
1675	CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1676	CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1677	CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1678
1679	! calculate lai and laimean (pixel including no_bio fraction)
1680	histvar(:)=zero
1681	DO ji = 1, kjpindex
1682	IF (SUM(veget_max(ji,:)) > zero) THEN
1683	DO jv=2,nvm
1684	histvar(ji) = histvar(ji) + veget_max(ji,jv)*lai(ji,jv)/ &
1685	SUM(veget_max(ji,:))
1686	END DO
1687	END IF
1688	END DO
1689	CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1690	CALL histwrite_p(hist_id, 'LAImean', kjit, histvar, kjpindex, index)
1691	CALL histwrite_p(hist_id, 'subli', kjit, vevapsno, kjpindex, index)
1692	CALL histwrite_p(hist_id, 'evapnu', kjit, vevapnu, kjpindex, index)
1693	CALL histwrite_p(hist_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
1694	CALL histwrite_p(hist_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
1695	CALL histwrite_p(hist_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
1696	CALL histwrite_p(hist_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
1697	CALL histwrite_p(hist_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
1698	CALL histwrite_p(hist_id, 'vbeta4', kjit, vbeta4, kjpindex, index)
1699	CALL histwrite_p(hist_id, 'vbeta5', kjit, vbeta5, kjpindex, index)
1700	CALL histwrite_p(hist_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
1701	CALL histwrite_p(hist_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
1702	CALL histwrite_p(hist_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
1703
1704	CALL histwrite_p(hist_id, 'snow', kjit, snow, kjpindex, index)
1705	CALL histwrite_p(hist_id, 'snowage', kjit, snow_age, kjpindex, index)
1706	CALL histwrite_p(hist_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
1707	CALL histwrite_p(hist_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
1708
1709	CALL histwrite_p(hist_id, 'grndflux', kjit, grndflux, kjpindex,index)
1710	CALL histwrite_p(hist_id, 'snowtemp',kjit,snowtemp,kjpindex*nsnow,indexsnow)
1711	CALL histwrite_p(hist_id, 'snowliq', kjit,snowliq,kjpindex*nsnow,indexsnow)
1712	CALL histwrite_p(hist_id, 'snowdz', kjit,snowdz,kjpindex*nsnow,indexsnow)
1713	CALL histwrite_p(hist_id, 'snowrho', kjit,snowrho,kjpindex*nsnow,indexsnow)
1714	CALL histwrite_p(hist_id, 'snowgrain',kjit,snowgrain,kjpindex*nsnow,indexsnow)
1715	CALL histwrite_p(hist_id, 'snowheat',kjit,snowheat,kjpindex*nsnow,indexsnow)
1716
1717	CALL histwrite_p(hist_id,'frac_snow_veg',kjit,frac_snow_veg,kjpindex,index)
1718	CALL histwrite_p(hist_id, 'frac_snow_nobio', kjit,frac_snow_nobio,kjpindex*nnobio, indexnobio)
1719	CALL histwrite_p(hist_id, 'pgflux',kjit,pgflux,kjpindex,index)
1720	CALL histwrite_p(hist_id, 'soiltile', kjit, soiltile, kjpindex*nstm, indexsoil)
1721
1722	CALL histwrite_p(hist_id, 'soilindex', kjit, REAL(njsc, r_std), kjpindex, index)
1723	CALL histwrite_p(hist_id, 'reinf_slope', kjit, reinf_slope, kjpindex, index)
1724	CALL histwrite_p(hist_id, 'k_litt', kjit, k_litt, kjpindex, index)
1725
1726	IF ( do_floodplains ) THEN
1727	CALL histwrite_p(hist_id, 'evapflo', kjit, vevapflo, kjpindex, index)
1728	CALL histwrite_p(hist_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
1729	ENDIF
1730
1731	CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)
1732	CALL histwrite_p(hist_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
1733	CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
1734
1735	IF ( ok_c13 ) THEN
1736	CALL histwrite_p(hist_id, 'delta_c13_assim', kjit, delta_c13_assim, &
1737	kjpindex*nvm, indexveg)
1738	CALL histwrite_p(hist_id, 'leaf_ci_out', kjit, leaf_ci_out, kjpindex*nvm, &
1739	indexveg)
1740	CALL histwrite_p(hist_id, 'gpp_day', kjit, gpp_day, kjpindex*nvm, &
1741	indexveg)
1742	ENDIF
1743
1744	histvar(:)=SUM(vevapwet(:,:),dim=2)
1745	CALL histwrite_p(hist_id, 'evspsblveg', kjit, histvar, kjpindex, index)
1746
1747	histvar(:)= vevapnu(:)+vevapsno(:)
1748	CALL histwrite_p(hist_id, 'evspsblsoi', kjit, histvar, kjpindex, index)
1749
1750	histvar(:)=SUM(transpir(:,:),dim=2)
1751	CALL histwrite_p(hist_id, 'tran', kjit, histvar, kjpindex, index)
1752
1753	histvar(:)= sum_treefrac(:)100contfrac(:)
1754	CALL histwrite_p(hist_id, 'treeFrac', kjit, histvar, kjpindex, index)
1755
1756	histvar(:)= (sum_grassfracC3(:)+sum_grassfracC4(:))100contfrac(:)
1757	CALL histwrite_p(hist_id, 'grassFrac', kjit, histvar, kjpindex, index)
1758
1759	histvar(:)= (sum_cropfracC3(:)+sum_cropfracC4(:))100contfrac(:)
1760	CALL histwrite_p(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)
1761
1762	histvar(:)=veget_max(:,1)100contfrac(:)
1763	CALL histwrite_p(hist_id, 'baresoilFrac', kjit, histvar, kjpindex, index)
1764
1765	histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)100contfrac(:)
1766	CALL histwrite_p(hist_id, 'residualFrac', kjit, histvar, kjpindex, index)
1767	ELSE
1768	! Write history file in ALMA format
1769	CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1770	CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1771	CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1772	CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1773	CALL histwrite_p(hist_id, 'ESoil', kjit, vevapnu, kjpindex, index)
1774	CALL histwrite_p(hist_id, 'EWater', kjit, vevapflo, kjpindex, index)
1775	CALL histwrite_p(hist_id, 'SWE', kjit, snow, kjpindex, index)
1776	histvar(:)=zero
1777	DO jv=1,nvm
1778	histvar(:) = histvar(:) + transpir(:,jv)
1779	ENDDO
1780	CALL histwrite_p(hist_id, 'TVeg', kjit, histvar, kjpindex, index)
1781	histvar(:)=zero
1782	DO jv=1,nvm
1783	histvar(:) = histvar(:) + vevapwet(:,jv)
1784	ENDDO
1785	CALL histwrite_p(hist_id, 'ECanop', kjit, histvar, kjpindex, index)
1786	CALL histwrite_p(hist_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
1787	!
1788	CALL histwrite_p(hist_id, 'Z0m', kjit, z0m, kjpindex, index)
1789	CALL histwrite_p(hist_id, 'Z0h', kjit, z0h, kjpindex, index)
1790	CALL histwrite_p(hist_id, 'EffectHeight', kjit, roughheight, kjpindex, index)
1791	!
1792	IF ( do_floodplains ) THEN
1793	CALL histwrite_p(hist_id, 'Qflood', kjit, vevapflo, kjpindex, index)
1794	CALL histwrite_p(hist_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
1795	ENDIF
1796
1797	CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)
1798	CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
1799	CALL histwrite_p(hist_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
1800
1801	ENDIF ! almaoutput
1802
1803	!! 13. Write additional output file with higher frequency
1804	IF ( hist2_id > 0 ) THEN
1805	IF ( .NOT. almaoutput ) THEN
1806	WRITE(numout,*)'vbeta, what is happening',vbeta
1807	! Write history file in IPSL-format
1808	CALL histwrite_p(hist2_id, 'tsol_rad', kjit, tsol_rad, kjpindex, index)
1809	CALL histwrite_p(hist2_id, 'qsurf', kjit, qsurf, kjpindex, index)
1810	CALL histwrite_p(hist2_id, 'albedo', kjit, albedo, kjpindex*2, indexalb)
1811	CALL histwrite_p(hist2_id, 'emis', kjit, emis, kjpindex, index)
1812	CALL histwrite_p(hist2_id, 'beta', kjit, vbeta, kjpindex, index)
1813	CALL histwrite_p(hist2_id, 'z0m', kjit, z0m, kjpindex, index)
1814	CALL histwrite_p(hist2_id, 'z0h', kjit, z0h, kjpindex, index)
1815	CALL histwrite_p(hist2_id, 'roughheight', kjit, roughheight, kjpindex, index)
1816	CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1817	CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1818	CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1819	CALL histwrite_p(hist2_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1820	CALL histwrite_p(hist2_id, 'subli', kjit, vevapsno, kjpindex, index)
1821	IF ( do_floodplains ) THEN
1822	CALL histwrite_p(hist2_id, 'vevapflo', kjit, vevapflo, kjpindex, index)
1823	CALL histwrite_p(hist2_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
1824	ENDIF
1825	CALL histwrite_p(hist2_id, 'vevapnu', kjit, vevapnu, kjpindex, index)
1826	CALL histwrite_p(hist2_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
1827	CALL histwrite_p(hist2_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
1828	CALL histwrite_p(hist2_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
1829	CALL histwrite_p(hist2_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
1830	CALL histwrite_p(hist2_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
1831	CALL histwrite_p(hist2_id, 'vbeta4', kjit, vbeta4, kjpindex, index)
1832	CALL histwrite_p(hist2_id, 'vbeta5', kjit, vbeta5, kjpindex, index)
1833	CALL histwrite_p(hist2_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
1834	CALL histwrite_p(hist2_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
1835	CALL histwrite_p(hist2_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
1836	CALL histwrite_p(hist2_id, 'snow', kjit, snow, kjpindex, index)
1837	CALL histwrite_p(hist2_id, 'snowage', kjit, snow_age, kjpindex, index)
1838	CALL histwrite_p(hist2_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
1839	CALL histwrite_p(hist2_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
1840
1841	IF (ok_hydrol_arch) THEN
1842	CALL histwrite_p(hist2_id, 'transpir_supply', kjit,transpir_supply, kjpindex*nvm, indexveg)
1843	ENDIF
1844
1845	CALL histwrite_p(hist2_id, 'soilindex', kjit, REAL(njsc, r_std), kjpindex, index)
1846	CALL histwrite_p(hist2_id, 'reinf_slope', kjit, reinf_slope, kjpindex, index)
1847
1848	CALL histwrite_p(hist2_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)
1849	CALL histwrite_p(hist2_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
1850	CALL histwrite_p(hist2_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
1851
1852	ELSE
1853	! Write history file in ALMA format
1854	CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1855	CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1856	CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1857	CALL histwrite_p(hist2_id, 'ESoil', kjit, vevapnu, kjpindex, index)
1858	IF ( do_floodplains ) THEN
1859	CALL histwrite_p(hist2_id, 'EWater', kjit, vevapflo, kjpindex, index)
1860	CALL histwrite_p(hist2_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
1861	ENDIF
1862	CALL histwrite_p(hist2_id, 'SWE', kjit, snow, kjpindex, index)
1863	histvar(:)=zero
1864	DO jv=1,nvm
1865	histvar(:) = histvar(:) + transpir(:,jv)
1866	ENDDO
1867	CALL histwrite_p(hist2_id, 'TVeg', kjit, histvar, kjpindex, index)
1868	histvar(:)=zero
1869	DO jv=1,nvm
1870	histvar(:) = histvar(:) + vevapwet(:,jv)
1871	ENDDO
1872	CALL histwrite_p(hist2_id, 'ECanop', kjit, histvar, kjpindex, index)
1873	CALL histwrite_p(hist2_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
1874	CALL histwrite_p(hist2_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
1875
1876	ENDIF ! almaoutput
1877	ENDIF ! hist2_id
1878
1879	END SUBROUTINE sechiba_main
1880
1881	!! =============================================================================================================================
1882	!! SUBROUTINE: sechiba_finalize
1883	!!
1884	!>\BRIEF Finalize all modules by calling their "_finalize" subroutines.
1885	!!
1886	!! DESCRIPTION: Finalize all modules by calling their "_finalize" subroutines. These subroutines will write variables to
1887	!! restart file.
1888	!!
1889	!! \n
1890	!_ ==============================================================================================================================
1891	!+++CHECK+++
1892	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
1893	! when running a larger domain. Needs to be corrected when implementing
1894	! a global use of the multi-layer energy budget.
1895	SUBROUTINE sechiba_finalize( &
1896	kjit, kjpij, kjpindex, index, rest_id, &
1897	tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad, &
1898	albedo &
1899	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean &
1900	)
1901	!+++++++++++
1902
1903	!! 0.1 Input variables
1904	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
1905	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
1906	!! (unitless)
1907	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
1908	!! (unitless)
1909	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
1910	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the pixels on the map.
1911	!! Sechiba uses a reduced grid excluding oceans
1912	!! ::index contains the indices of the
1913	!! terrestrial pixels only! (unitless)
1914	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tsol_rad !! Radiative surface temperature
1915	!! @tex $(W m^{-2})$ @endtex
1916	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapp !! Total of evaporation
1917	!! @tex $(kg m^{-2} days^{-1})$ @endtex
1918	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxsens !! Sensible heat flux
1919	!! @tex $(W m^{-2})$ @endtex
1920	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxlat !! Latent heat flux
1921	!! @tex $(W m^{-2})$ @endtex
1922	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tq_cdrag !! Surface drag coefficient (-)
1923	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: albedo !! Surface albedo for visible and near-infrared (unitless, 0-1)
1924	!+++CHECK+++
1925	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
1926	! when running a larger domain. Needs to be corrected when implementing
1927	! a global use of the multi-layer energy budget.
1928	!!$ REAL(r_std),DIMENSION(nlevels_tot) :: Light_Abs_Tot_mean !! total light absorption for a given canopy level
1929	!!$ REAL(r_std),DIMENSION(nlevels_tot) :: Light_Alb_Tot_mean !! total albedo for a given level
1930	!+++++++++++
1931
1932	!! 0.2 Local variables
1933	INTEGER(i_std) :: ji, jv !! Index (unitless)
1934	REAL(r_std), DIMENSION(kjpindex) :: histvar !! Computations for history files (unitless)
1935	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
1936
1937	! ==============================================================================================================================
1938
1939	!! Write restart file for the next simulation from SECHIBA and other modules
1940	IF (printlev_loc>=3) WRITE (numout,*) 'Start sechiba_finalize for writing restart files'
1941
1942	!! 1. Call diffuco_finalize to write restart files
1943	CALL diffuco_finalize (kjit, kjpindex, rest_id, rstruct )
1944
1945
1946	!! 2. Call energy budget to write restart files
1947	IF(ok_mleb) THEN
1948	! Multi-layer energy budget
1949	CALL mleb_finalize (kjit, kjpindex, rest_id, &
1950	evapot, evapot_corr, temp_sol, tsol_rad, &
1951	qsurf, fluxsens, fluxlat, vevapp, &
1952	u_speed, z_array_out, &
1953	max_height_store )
1954	ELSE
1955	! Energy budget using enerbil module
1956	CALL enerbil_finalize(kjit, kjpindex, rest_id, &
1957	evapot, evapot_corr, temp_sol, tsol_rad, &
1958	qsurf, fluxsens, fluxlat, vevapp )
1959	ENDIF
1960
1961
1962	!! 3. Call hydrology to write restart files
1963	CALL hydrol_finalize( kjit, kjpindex, rest_id, vegstress, &
1964	qsintveg, humrel, snow, snow_age, snow_nobio, &
1965	snow_nobio_age, snowrho, snowtemp, snowdz, &
1966	snowheat, snowgrain, drysoil_frac, evap_bare_lim, &
1967	evap_bare_lim_ns, swc, ksoil, root_profile, &
1968	us)
1969
1970
1971	!! 4. Call condveg to write surface variables to restart files
1972	!+++CHECK+++
1973	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
1974	! when running a larger domain. Needs to be corrected when implementing
1975	! a global use of the multi-layer energy budget.
1976	CALL condveg_finalize (kjit, kjpindex, rest_id, &
1977	z0m, z0h, roughheight, &
1978	albedo, Light_Abs_Tot, Light_Tran_Tot, &
1979	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
1980	laieff_isotrop)
1981	!+++++++++++
1982
1983	!! 5. Call soil thermodynamic to write restart files
1984	CALL thermosoil_finalize (kjit, kjpindex, rest_id, gtemp, &
1985	soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
1986
1987	!! 6. Add river routing to restart files
1988	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
1989	!! 6.1 Call river routing to write restart files
1990	CALL routing_finalize( kjit, kjpindex, rest_id, flood_frac, flood_res )
1991	ELSE
1992	!! 6.2 No routing, set variables to zero
1993	reinfiltration(:) = zero
1994	returnflow(:) = zero
1995	irrigation(:) = zero
1996	flood_frac(:) = zero
1997	flood_res(:) = zero
1998	ENDIF
1999
2000	!! 7. Call slowproc_main to add 'daily' and annual variables to restart file
2001	CALL slowproc_finalize (kjit, kjpindex, rest_id, index, &
2002	njsc, veget, frac_nobio, &
2003	veget_max, reinf_slope, ks, nvan, &
2004	avan, mcr, mcs, mcfc, &
2005	mcw, assim_param, frac_age, &
2006	heat_Zimov,altmax, depth_organic_soil, &
2007	circ_class_biomass, circ_class_n, &
2008	lai_per_level, laieff_fit, loss_gain, &
2009	veget_max_new, frac_nobio_new)
2010
2011	!! 8. Write variables from sechiba module to restart file
2012	!! Here are only the variables set which are not already restarted from other modules.
2013	CALL restput_p (rest_id, 'coszang_noon', nbp_glo, 1, 1, kjit, coszang_noon, 'scatter', nbp_glo, index_g)
2014
2015	IF (printlev_loc>=3) WRITE (numout,*) 'sechiba_finalize done'
2016
2017	END SUBROUTINE sechiba_finalize
2018
2019
2020	!! ==============================================================================================================================\n
2021	!! SUBROUTINE : sechiba_init
2022	!!
2023	!>\BRIEF Dynamic allocation of the variables, the dimensions of the
2024	!! variables are determined by user-specified settings.
2025	!!
2026	!! DESCRIPTION : The domain size (:: kjpindex) is used to allocate the correct
2027	!! dimensions to all variables in sechiba. Depending on the variable, its
2028	!! dimensions are also determined by the number of PFT's (::nvm), number of
2029	!! soil types (::nstm), number of non-vegetative surface types (::nnobio),
2030	!! number of soil levels (::ngrnd), number of soil layers in the hydrological
2031	!! model (i.e. cwrr) (::nslm). Values for these variables are set in
2032	!! constantes_soil.f90 and constantes_veg.f90.\n
2033	!!
2034	!! Memory is allocated for all Sechiba variables and new indexing tables
2035	!! are build making use of both (::kjpij) and (::kjpindex). New indexing tables
2036	!! are needed because a single pixel can contain several PFTs, soil types, etc.
2037	!! The new indexing tables have separate indices for the different
2038	!! PFTs, soil types, etc.\n
2039	!!
2040	!! RECENT CHANGE(S): None
2041	!!
2042	!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
2043	!! variables. However, the routine allocates memory and builds new indexing
2044	!! variables for later use.
2045	!!
2046	!! REFERENCE(S) : None
2047	!!
2048	!! FLOWCHART : None
2049	!! \n
2050	!_ ================================================================================================================================
2051
2052	SUBROUTINE sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)
2053
2054	!! 0.1 Input variables
2055
2056	INTEGER(i_std), INTENT (in) :: kjit !! Time step number (unitless)
2057	INTEGER(i_std), INTENT (in) :: kjpij !! Total size of the un-compressed grid (unitless)
2058	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
2059	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
2060	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless)
2061	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates (latitude,longitude)
2062	!! for pixels (degrees)
2063	!! 0.2 Output variables
2064
2065	!! 0.3 Modified variables
2066
2067	!! 0.4 Local variables
2068
2069	INTEGER(i_std) :: ier !! Check errors in memory allocation (unitless)
2070	INTEGER(i_std) :: ji,jv,ilevel,ipts !! Indeces (unitless)
2071	CHARACTER(LEN=80) :: var_name !! To store variables names in restart file
2072	!_ ==============================================================================================================================
2073
2074	!! 1. Initialize variables
2075
2076	! Debug
2077	! It is good to leave this in here. It is only written out once,
2078	! it takes almost no time, and it's necessary for identifying a
2079	! problem pixel.
2080	DO ipts=1,kjpindex
2081	WRITE(numout,'(A,I6,10F20.10)') 'pixel number to lat/lon: ',ipts,lalo(ipts,1:2)
2082	ENDDO
2083	!-
2084
2085	! Dynamic allocation with user-specified dimensions on first call
2086	IF (l_first_sechiba) THEN
2087	l_first_sechiba=.FALSE.
2088	ELSE
2089	CALL ipslerr_p(3,'sechiba_init',' l_first_sechiba false . we stop ','','')
2090	ENDIF
2091
2092	!! 1.1 Initialize 3D vegetation indexation table
2093	ALLOCATE (indexveg(kjpindex*nvm),stat=ier)
2094	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexveg','','')
2095
2096	ALLOCATE (indexlai(kjpindex*(nlai+1)),stat=ier)
2097	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlai','','')
2098
2099	ALLOCATE (indexsoil(kjpindex*nstm),stat=ier)
2100	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsoil','','')
2101
2102	ALLOCATE (indexnobio(kjpindex*nnobio),stat=ier)
2103	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnobio','','')
2104
2105	ALLOCATE (indexgrnd(kjpindex*ngrnd),stat=ier)
2106	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexgrnd','','')
2107
2108	ALLOCATE (indexsnow(kjpindex*nsnow),stat=ier)
2109	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsnow','','')
2110
2111	ALLOCATE (indexlayer(kjpindex*nslm),stat=ier)
2112	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlayer','','')
2113
2114	ALLOCATE (indexnslm(kjpindex*nslm),stat=ier)
2115	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnslm','','')
2116
2117	ALLOCATE (indexalb(kjpindex*2),stat=ier)
2118	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexalb','','')
2119
2120	!! 1.2 Initialize 1D array allocation with restartable value
2121	ALLOCATE (flood_res(kjpindex),stat=ier)
2122	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for flood_res','','')
2123	flood_res(:) = undef_sechiba
2124
2125	ALLOCATE (flood_frac(kjpindex),stat=ier)
2126	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for kjpindex','','')
2127	flood_frac(:) = undef_sechiba
2128
2129	ALLOCATE (snow(kjpindex),stat=ier)
2130	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow','','')
2131	snow(:) = undef_sechiba
2132
2133	ALLOCATE (snow_age(kjpindex),stat=ier)
2134	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_age','','')
2135	snow_age(:) = undef_sechiba
2136
2137	ALLOCATE (drysoil_frac(kjpindex),stat=ier)
2138	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drysoil_frac','','')
2139
2140	ALLOCATE (evap_bare_lim(kjpindex),stat=ier)
2141	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim','','')
2142
2143	ALLOCATE (evap_bare_lim_ns(kjpindex,nstm),stat=ier)
2144	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim_ns','','')
2145
2146	ALLOCATE (evapot(kjpindex),stat=ier)
2147	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot','','')
2148	evapot(:) = undef_sechiba
2149
2150	ALLOCATE (evapot_corr(kjpindex),stat=ier)
2151	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot_corr','','')
2152
2153	ALLOCATE (humrel(kjpindex,nvm),stat=ier)
2154	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for humrel','','')
2155	humrel(:,:) = undef_sechiba
2156
2157	ALLOCATE (vegstress(kjpindex,nvm),stat=ier)
2158	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vegstress','','')
2159	vegstress(:,:) = undef_sechiba
2160
2161	ALLOCATE (njsc(kjpindex),stat=ier)
2162	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for njsc','','')
2163	njsc(:)= undef_int
2164
2165	ALLOCATE (soiltile(kjpindex,nstm),stat=ier)
2166	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soiltile','','')
2167
2168	ALLOCATE (fraclut(kjpindex,nlut),stat=ier)
2169	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fraclut','','')
2170
2171	ALLOCATE (nwdFraclut(kjpindex,nlut),stat=ier)
2172	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for nwdFraclut','','')
2173
2174	ALLOCATE (reinf_slope(kjpindex),stat=ier)
2175	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope','','')
2176
2177	ALLOCATE (ks(kjpindex),stat=ier)
2178	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ks','','')
2179
2180	ALLOCATE (nvan(kjpindex),stat=ier)
2181	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for nvan ','','')
2182
2183	ALLOCATE (avan(kjpindex),stat=ier)
2184	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for avan','','')
2185
2186	ALLOCATE (mcr(kjpindex),stat=ier)
2187	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcr','','')
2188
2189	ALLOCATE (mcs(kjpindex),stat=ier)
2190	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcs','','')
2191
2192	ALLOCATE (mcfc(kjpindex),stat=ier)
2193	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcfc','','')
2194
2195	ALLOCATE (mcw(kjpindex),stat=ier)
2196	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcw','','')
2197
2198	ALLOCATE (vbeta1(kjpindex),stat=ier)
2199	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta1','','')
2200
2201	ALLOCATE (us(kjpindex,nvm,nstm,nslm),stat=ier)
2202	IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable us','','')
2203
2204	ALLOCATE (vbeta4(kjpindex),stat=ier)
2205	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta4','','')
2206
2207	ALLOCATE (vbeta5(kjpindex),stat=ier)
2208	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta5','','')
2209
2210	ALLOCATE (soilcap(kjpindex),stat=ier)
2211	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilcap','','')
2212
2213	ALLOCATE (soilflx(kjpindex),stat=ier)
2214	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilflx','','')
2215
2216	ALLOCATE (temp_sol(kjpindex),stat=ier)
2217	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol','','')
2218	temp_sol(:) = undef_sechiba
2219
2220	ALLOCATE (qsurf(kjpindex),stat=ier)
2221	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsurf','','')
2222	qsurf(:) = undef_sechiba
2223
2224	!! 1.3 Initialize 2D array allocation with restartable value
2225	ALLOCATE (qsintveg(kjpindex,nvm),stat=ier)
2226	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintveg','','')
2227	qsintveg(:,:) = undef_sechiba
2228
2229	ALLOCATE (vbeta2(kjpindex,nvm),stat=ier)
2230	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta2','','')
2231
2232	ALLOCATE (vbeta3(kjpindex,nvm),stat=ier)
2233	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3','','')
2234
2235	ALLOCATE (vbeta3pot(kjpindex,nvm),stat=ier)
2236	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3pot','','')
2237
2238	ALLOCATE (gsmean(kjpindex,nvm),stat=ier)
2239	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gsmean','','')
2240
2241	ALLOCATE (cimean(kjpindex,nvm),stat=ier)
2242	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cimean','','')
2243
2244	ALLOCATE (gpp(kjpindex,nvm),stat=ier)
2245	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gpp','','')
2246	gpp(:,:) = undef_sechiba
2247
2248	ALLOCATE (temp_growth(kjpindex),stat=ier)
2249	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_growth','','')
2250	temp_growth(:) = undef_sechiba
2251
2252	ALLOCATE (veget(kjpindex,nvm),stat=ier)
2253	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget','','')
2254	veget(:,:)=undef_sechiba
2255
2256	ALLOCATE (veget_max(kjpindex,nvm),stat=ier)
2257	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget_max','','')
2258
2259	ALLOCATE (tot_bare_soil(kjpindex),stat=ier)
2260	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_bare_soil','','')
2261
2262	ALLOCATE (lai(kjpindex,nvm),stat=ier)
2263	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lai','','')
2264	lai(:,:)=undef_sechiba
2265
2266	ALLOCATE (laieff_fit(kjpindex,nvm,nlevels_tot),stat=ier)
2267	IF (ier.NE.0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for laieff_fit','','')
2268	CALL laieff_type_init(kjpindex, nlevels_tot, laieff_fit)
2269
2270	ALLOCATE (frac_age(kjpindex,nvm,nleafages),stat=ier)
2271	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_age','','')
2272	frac_age(:,:,:)=undef_sechiba
2273
2274	ALLOCATE (height(kjpindex,nvm),stat=ier)
2275	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for height','','')
2276	height(:,:)=undef_sechiba
2277
2278	ALLOCATE (height_dom(kjpindex,nvm),stat=ier)
2279	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for height_dom','','')
2280	height_dom(:,:)=undef_sechiba
2281
2282	ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
2283	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_nobio','','')
2284	frac_nobio(:,:) = undef_sechiba
2285
2286	ALLOCATE (snow_nobio(kjpindex,nnobio),stat=ier)
2287	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio','','')
2288	snow_nobio(:,:) = undef_sechiba
2289
2290	ALLOCATE (snow_nobio_age(kjpindex,nnobio),stat=ier)
2291	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio_age','','')
2292	snow_nobio_age(:,:) = undef_sechiba
2293
2294	ALLOCATE (assim_param(kjpindex,nvm,npco2),stat=ier)
2295	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for assim_param','','')
2296
2297	!! 1.4 Initialize 1D array allocation
2298	ALLOCATE (vevapflo(kjpindex),stat=ier)
2299	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapflo','','')
2300	vevapflo(:)=zero
2301
2302	ALLOCATE (vevapsno(kjpindex),stat=ier)
2303	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapsno','','')
2304
2305	ALLOCATE (vevapnu(kjpindex),stat=ier)
2306	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapnu','','')
2307
2308	ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
2309	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for totfrac_nobio','','')
2310
2311	ALLOCATE (floodout(kjpindex),stat=ier)
2312	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for floodout','','')
2313
2314	ALLOCATE (runoff(kjpindex),stat=ier)
2315	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for runoff','','')
2316
2317	ALLOCATE (drainage(kjpindex),stat=ier)
2318	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drainage','','')
2319
2320	ALLOCATE (returnflow(kjpindex),stat=ier)
2321	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for returnflow','','')
2322	returnflow(:) = zero
2323
2324	ALLOCATE (reinfiltration(kjpindex),stat=ier)
2325	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinfiltration','','')
2326	reinfiltration(:) = zero
2327
2328	ALLOCATE (irrigation(kjpindex),stat=ier)
2329	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for irrigation','','')
2330	irrigation(:) = zero
2331
2332	ALLOCATE (z0h(kjpindex),stat=ier)
2333	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0h','','')
2334
2335	ALLOCATE (z0m(kjpindex),stat=ier)
2336	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0m','','')
2337
2338	ALLOCATE (roughheight(kjpindex),stat=ier)
2339	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for roughheight','','')
2340
2341	ALLOCATE (emis(kjpindex),stat=ier)
2342	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for emis','','')
2343
2344	ALLOCATE (tot_melt(kjpindex),stat=ier)
2345	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_melt','','')
2346
2347	ALLOCATE (vbeta(kjpindex),stat=ier)
2348	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta','','')
2349
2350	ALLOCATE (rau(kjpindex),stat=ier)
2351	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rau','','')
2352
2353	ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
2354	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for deadleaf_cover','','')
2355
2356	ALLOCATE (stempdiag(kjpindex, nslm),stat=ier)
2357	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for stempdiag','','')
2358
2359	ALLOCATE (co2_to_bm(kjpindex,nvm),stat=ier)
2360	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for co2_to_bm','','')
2361
2362	ALLOCATE (shumdiag(kjpindex,nslm),stat=ier)
2363	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag','','')
2364
2365	ALLOCATE (shumdiag_perma(kjpindex,nslm),stat=ier)
2366	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag_perma','','')
2367
2368	ALLOCATE (litterhumdiag(kjpindex),stat=ier)
2369	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for litterhumdiag','','')
2370
2371	ALLOCATE (ptnlev1(kjpindex),stat=ier)
2372	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ptnlev1','','')
2373
2374	ALLOCATE (k_litt(kjpindex),stat=ier)
2375	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for k_litt','','')
2376
2377	ALLOCATE (netrad(kjpindex),stat=ier)
2378	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for netrad','','')
2379
2380	ALLOCATE (lwabs(kjpindex),stat=ier)
2381	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lwabs','','')
2382
2383	ALLOCATE (lwnet(kjpindex),stat=ier)
2384	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lwnet','','')
2385
2386	ALLOCATE (fluxsubli(kjpindex),stat=ier)
2387	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fluxsubli','','')
2388
2389
2390	!! 1.5 Initialize 2D array allocation
2391	ALLOCATE (vevapwet(kjpindex,nvm),stat=ier)
2392	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapwet','','')
2393	vevapwet(:,:)=undef_sechiba
2394
2395	ALLOCATE (transpir(kjpindex,nvm),stat=ier)
2396	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpir','','')
2397
2398	ALLOCATE (transpot(kjpindex,nvm),stat=ier)
2399	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpot','','')
2400
2401	ALLOCATE (transpir_mod(kjpindex,nvm),stat=ier)
2402	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init','Pb in alloc for transpir_mod','','')
2403	transpir_mod(:,:) = zero
2404
2405
2406	ALLOCATE (stressed(kjpindex,nvm),stat=ier)
2407	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init','Pb in alloc for stressed','','')
2408	stressed(:,:) = zero
2409
2410	ALLOCATE (unstressed(kjpindex,nvm),stat=ier)
2411	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init','Pb in alloc for unstressed','','')
2412	unstressed(:,:) = zero
2413
2414	ALLOCATE (transpir_supply(kjpindex,nvm),stat=ier)
2415	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init','Pb in alloc for transpir_supply','','')
2416	transpir_supply(:,:) = zero
2417
2418	ALLOCATE (vir_transpir_supply(kjpindex,nvm),stat=ier)
2419	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init','Pb in alloc for vir_transpir_supply','','')
2420	vir_transpir_supply(:,:) = zero
2421
2422	ALLOCATE (transpir_supply_column(nlevels_tot,kjpindex,nvm),stat=ier)
2423	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init','Pb in alloc for transpir_supply_column','','')
2424	transpir_supply_column(:,:,:) = zero
2425
2426	ALLOCATE (vessel_loss(kjpindex,nvm),stat=ier)
2427	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init','Pb in alloc for vessel_loss','','')
2428
2429	ALLOCATE (e_frac(kjpindex,nvm,nslm,nstm),stat=ier)
2430	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init','Pb in alloc for e_frac','','')
2431	e_frac(:,:,:,:) = zero
2432
2433	ALLOCATE (qsintmax(kjpindex,nvm),stat=ier)
2434	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintmax','','')
2435
2436	ALLOCATE (rveget(kjpindex,nvm),stat=ier)
2437	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rveget','','')
2438
2439	ALLOCATE (rstruct(kjpindex,nvm),stat=ier)
2440	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rstruct','','')
2441
2442	ALLOCATE (pgflux(kjpindex),stat=ier)
2443	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for pgflux','','')
2444	pgflux(:)= 0.0
2445
2446	ALLOCATE (cgrnd_snow(kjpindex,nsnow),stat=ier)
2447	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cgrnd_snow','','')
2448	cgrnd_snow(:,:) = 0
2449
2450	ALLOCATE (dgrnd_snow(kjpindex,nsnow),stat=ier)
2451	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for dgrnd_snow','','')
2452	dgrnd_snow(:,:) = 0
2453
2454	ALLOCATE (lambda_snow(kjpindex),stat=ier)
2455	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lambda_snow','','')
2456	lambda_snow(:) = 0
2457
2458	ALLOCATE (temp_sol_add(kjpindex),stat=ier)
2459	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol_add','','')
2460
2461	ALLOCATE (qsol_sat_new(kjpindex),stat=ier)
2462	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsol_sat_new','','')
2463
2464	ALLOCATE (qair_new(kjpindex),stat=ier)
2465	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qair_new','','')
2466
2467	ALLOCATE (gtemp(kjpindex),stat=ier)
2468	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gtemp','','')
2469
2470	ALLOCATE (frac_snow_veg(kjpindex),stat=ier)
2471	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_snow_veg','','')
2472
2473	ALLOCATE (frac_snow_nobio(kjpindex,nnobio),stat=ier)
2474	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_snow_nobio','','')
2475
2476	ALLOCATE (snowrho(kjpindex,nsnow),stat=ier)
2477	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowrho','','')
2478
2479	ALLOCATE (snowheat(kjpindex,nsnow),stat=ier)
2480	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowheat','','')
2481
2482	ALLOCATE (snowgrain(kjpindex,nsnow),stat=ier)
2483	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowgrain','','')
2484
2485	ALLOCATE (snowtemp(kjpindex,nsnow),stat=ier)
2486	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowtemp','','')
2487
2488	ALLOCATE (snowdz(kjpindex,nsnow),stat=ier)
2489	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowdz','','')
2490
2491	ALLOCATE(max_height_store(kjpindex,nvm),stat=ier)
2492	IF (ier.NE.0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for max_height_store','','')
2493
2494	ALLOCATE (warnings(kjpindex,nvm,nwarns),stat=ier)
2495	IF (ier.NE.0) THEN
2496	WRITE (numout,*) ' error in warnings allocation. We stop. We need kjpindex x nvm x nwarns words = ',&
2497	& kjpindex,' x ' ,nvm, ' x ',nwarns,' = ',kjpindexnvmnwarns
2498	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2499	END IF
2500	! We aren't restarting the warnings at the moment.
2501	warnings(:,:,:)=zero
2502
2503	!! 1.5b Initialize 3D array allocation for albedo
2504	ALLOCATE (Light_Abs_Tot(kjpindex,nvm,nlevels_tot),stat=ier)
2505	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init', &
2506	'Pb in alloc for Light_Abs_Tot','','')
2507
2508	ALLOCATE (Light_Tran_Tot(kjpindex,nvm,nlevels_tot),stat=ier)
2509	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init', &
2510	'Pb in alloc for Light_Tran_Tot','','')
2511
2512	ALLOCATE (laieff_isotrop(kjpindex,nlevels_tot,nvm),stat=ier)
2513	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init', &
2514	'Pb in alloc for laieff_isotrop','','')
2515
2516	ALLOCATE (mc_layh(kjpindex, nslm),stat=ier)
2517	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mc_layh','','')
2518
2519	ALLOCATE (mcl_layh(kjpindex, nslm),stat=ier)
2520	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcl_layh','','')
2521
2522	ALLOCATE (soilmoist(kjpindex, nslm),stat=ier)
2523	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilmoist','','')
2524
2525	ALLOCATE (mcl_layh_s(kjpindex, nslm, nstm),stat=ier)
2526	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcl_layh_s','','')
2527
2528	ALLOCATE (mc_layh_s(kjpindex, nslm, nstm),stat=ier)
2529	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mc_layh_s','','')
2530
2531	ALLOCATE (soilmoist_s(kjpindex, nslm, nstm),stat=ier)
2532	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilmoist_s','','')
2533
2534	ALLOCATE(tdeep(kjpindex,ngrnd,nvm),stat=ier)
2535	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tdeep','','')
2536
2537	ALLOCATE(hsdeep(kjpindex,ngrnd,nvm),stat=ier)
2538	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for hsdeep','','')
2539
2540	ALLOCATE(heat_Zimov(kjpindex,ngrnd,nvm),stat=ier)
2541	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for heat_Zimov','','')
2542
2543	ALLOCATE(sfluxCH4_deep(kjpindex),stat=ier)
2544	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for sfluxCH4_deep','','')
2545
2546	ALLOCATE(sfluxCO2_deep(kjpindex),stat=ier)
2547	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for sfluxCO2_deep','','')
2548
2549	ALLOCATE(som_total(kjpindex,ngrnd,nvm,nelements),stat=ier)
2550	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for som_total','','')
2551
2552	ALLOCATE (altmax(kjpindex,nvm),stat=ier)
2553	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for altmax','','')
2554	ALLOCATE(depth_organic_soil(kjpindex),stat=ier)
2555	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for depth_organic_soil','','')
2556
2557	!! 1.6 Initialize indexing table for the vegetation fields.
2558	! In SECHIBA we work on reduced grids but to store in the full 3D filed vegetation variable
2559	! we need another index table : indexveg, indexsoil, indexnobio and indexgrnd
2560	DO ji = 1, kjpindex
2561	!
2562	DO jv = 1, nlai+1
2563	indexlai((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
2564	ENDDO
2565	!
2566	DO jv = 1, nvm
2567	indexveg((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
2568	ENDDO
2569	!
2570	DO jv = 1, nstm
2571	indexsoil((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
2572	ENDDO
2573	!
2574	DO jv = 1, nnobio
2575	indexnobio((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
2576	ENDDO
2577	!
2578	DO jv = 1, ngrnd
2579	indexgrnd((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
2580	ENDDO
2581	!
2582	DO jv = 1, nsnow
2583	indexsnow((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
2584	ENDDO
2585
2586	DO jv = 1, nslm
2587	indexnslm((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
2588	ENDDO
2589
2590	DO jv = 1, nslm
2591	indexlayer((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
2592	ENDDO
2593	!
2594	DO jv = 1, 2
2595	indexalb((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
2596	ENDDO
2597	!
2598	ENDDO
2599
2600	! 1.7 now we create and initialize some plant variables that need to be passed
2601	! around sechiba and stomate (from the merge, there might be more to
2602	! put here)
2603	ALLOCATE(nlevels_loc(kjpindex),stat=ier)
2604	IF (ier .NE. 0) CALL ipslerr_p (3,'sechiba_init','Pb in alloc for nlevels_loc','','')
2605	nlevels_loc(:) = val_exp
2606
2607	ALLOCATE (z_array_out(kjpindex,nvm,ncirc,nlevels_tot),STAT=ier)
2608	IF (ier.NE.0) CALL ipslerr_p (3,'sechiba_init','Pb in alloc for z_array_out','','')
2609
2610	ALLOCATE (profile_vbeta3(kjpindex,nvm,nlevels_tot),STAT=ier)
2611	IF (ier.NE.0) THEN
2612	WRITE (numout,*) ' error in profile_vbeta3 allocation.'
2613	WRITE (numout,) 'We stop. We need kjpindexnvm words =',kjpindexnvmnlevels_tot
2614	STOP
2615	END IF
2616
2617	ALLOCATE (profile_rveget(kjpindex,nvm,nlevels_tot),STAT=ier)
2618	IF (ier.NE.0) THEN
2619	WRITE (numout,*) ' error in profile_rveget allocation.'
2620	WRITE (numout,) 'We stop. We need kjpindexnvm words =',kjpindexnvmnlevels_tot
2621	STOP
2622	END IF
2623
2624	ALLOCATE (delta_c13_assim(kjpindex,nvm),stat=ier)
2625	IF (ier.NE.0) THEN
2626	WRITE (numout,*) ' error in delta_c13_assim allocation. We stop. We need kjpindex x nvm words = ',&
2627	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
2628	STOP 'sechiba_init'
2629	END IF
2630	delta_c13_assim(:,:)=undef_sechiba
2631
2632	ALLOCATE (leaf_ci_out(kjpindex,nvm),stat=ier)
2633	IF (ier.NE.0) THEN
2634	WRITE (numout,*) ' error in leaf_ci_out allocation. We stop. We need kjpindex x nvm words = ',&
2635	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
2636	STOP 'sechiba_init'
2637	END IF
2638	leaf_ci_out(:,:)=undef_sechiba
2639
2640	ALLOCATE (gpp_day(kjpindex,nvm),stat=ier)
2641	IF (ier.NE.0) THEN
2642	WRITE (numout,*) ' error in gpp_day allocation. We stop. We need kjpindex x nvm words = ',&
2643	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
2644	STOP 'sechiba_init'
2645	END IF
2646	gpp_day(:,:)=undef_sechiba
2647
2648	ALLOCATE(circ_class_n(kjpindex,nvm,ncirc),stat=ier)
2649	IF (ier .NE. 0) THEN
2650	WRITE(numout,) 'Memory allocation error for circ_class_n. We stop. We need kjpindexnvm*ncirc words', &
2651	& kjpindex,nvm,ncirc
2652	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2653	ENDIF
2654	circ_class_n(:,:,:) = val_exp
2655
2656	ALLOCATE(circ_class_biomass(kjpindex,nvm,ncirc,nparts,nelements),stat=ier)
2657	IF (ier .NE. 0) THEN
2658	WRITE(numout,*) 'Memory allocation error for circ_class_biomass.'
2659	WRITE(numout,) 'We stop. We need kjpindexnvmnpartsncirc*nelmements words', &
2660	& kjpindex,nvm,ncirc,nparts,nelements
2661	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2662	ENDIF
2663	circ_class_biomass(:,:,:,:,:) = val_exp
2664
2665	ALLOCATE(lai_per_level(kjpindex,nvm,nlevels_tot),stat=ier)
2666	IF (ier .NE. 0) THEN
2667	WRITE(numout,*) 'Memory allocation error for lai_per_level. We stop. '//&
2668	'We need kjpindexnvmnlevels_tot words', &
2669	& kjpindex,nvm,nlevels_tot
2670	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2671	ENDIF
2672	lai_per_level(:,:,:)=val_exp
2673
2674	ALLOCATE (loss_gain(kjpindex,nvm),stat=ier)
2675	IF (ier.NE.0) THEN
2676	WRITE (numout,*) ' error in loss_gain allocation. We stop. We need kjpindex x nvm words = ',&
2677	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
2678	STOP 'sechiba_init'
2679	END IF
2680	loss_gain(:,:)=undef_sechiba
2681
2682	ALLOCATE (veget_max_new(kjpindex,nvm),stat=ier)
2683	IF (ier.NE.0) THEN
2684	WRITE (numout,*) ' error in veget_max_new allocation. We stop. We need kjpindex x nvm words = ',&
2685	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
2686	STOP 'sechiba_init'
2687	END IF
2688	veget_max_new(:,:)=undef_sechiba
2689
2690	ALLOCATE (frac_nobio_new(kjpindex,nnobio),stat=ier)
2691	IF (ier.NE.0) THEN
2692	WRITE (numout,*) ' error in frac_nobio_new allocation. We stop. We need kjpindex x nnobio words = ',&
2693	& kjpindex,' x ' ,nnobio, ' = ',kjpindex*nnobio
2694	STOP 'sechiba_init'
2695	END IF
2696	frac_nobio_new(:,:)=undef_sechiba
2697
2698	! multilayer allocations
2699	ALLOCATE (u_speed(jnlvls),stat=ier)
2700	IF (ier.NE.0) THEN
2701	WRITE (numout,*) ' error in u_speed allocation. We stop. We need kjpindex words = ',jnlvls
2702	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2703	END IF
2704
2705	ALLOCATE (swc(kjpindex,nslm,nstm))
2706	IF (ier.NE.0) THEN
2707	WRITE (numout,) ' error in swc allocation. We stop. We need kjpindex words = ',kjpindexnslm*nstm
2708	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2709	END IF
2710
2711	ALLOCATE (ksoil(kjpindex,nslm,nstm))
2712	IF (ier.NE.0) THEN
2713	WRITE (numout,) ' error in ksoil allocation. We stop. We need kjpindexwords = ',kjpindexnslm*nstm
2714	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2715	END IF
2716
2717	ALLOCATE (root_profile(kjpindex,nvm,nslm,nroot_prof))
2718	IF (ier.NE.0) THEN
2719	WRITE (numout,) ' error in root_profile allocation. We stop. We need kjpindexwords = ',kjpindexnvmnslmnroot_prof
2720	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2721	END IF
2722
2723	ALLOCATE (root_depth(kjpindex,nvm,ndepths))
2724	IF (ier.NE.0) THEN
2725	WRITE (numout,) ' error in root_depth allocation. We stop. We need kjpindexwords = ',kjpindexnvm*ndepths
2726	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2727	END IF
2728
2729	ALLOCATE (coszang_noon(kjpindex))
2730	IF (ier.NE.0) THEN
2731	WRITE (numout,*) ' error in coszang_noon allocation. We stop. We need kjpindexwords = ',kjpindex
2732	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2733	END IF
2734
2735	! Read coszang_noon from restart file
2736	var_name= 'coszang_noon'
2737	CALL ioconf_setatt_p('UNITS', '-')
2738	CALL ioconf_setatt_p('LONG_NAME','Solar angel at noon')
2739	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, &
2740	.TRUE., coszang_noon, "gather", nbp_glo, index_g)
2741	IF ( ALL( coszang_noon(:) .EQ. val_exp ) ) THEN
2742	! coszang_noon was not found in restart file.
2743	WRITE (numout,*) 'coszang_noon was not found in restart file. Call solarang_noon to calculate it.'
2744	CALL solarang_noon(kjpindex, lalo, coszang_noon)
2745	ELSEIF ( ANY( ABS(coszang_noon(:)) .GT. 2 ) ) THEN
2746	! This is a test if coszang_noon is out of range. This can be the case if
2747	! there is a mismatch in the land-sea mask. coszang_noon is the first
2748	! restart variable read in ORCHIDEE.
2749	CALL ipslerr_p(3, 'sechiba_init', 'coszang_noon from restart file is out of range.',&
2750	'There is probably a mismatch in the land/sea mask in the model and in the restart file.',&
2751	'Start the model again without restart files for ORCHIDEE.')
2752	END IF
2753
2754	ALLOCATE (Pgap_cumul(kjpindex,nvm,nlevels_tot))
2755	IF (ier.NE.0) THEN
2756	WRITE (numout,*) ' error in Pgap_cumul allocation. We stop. We need kjpindexwords = ', kjpindex,nvm,nlevels_tot
2757	CALL ipslerr_p (3,'sechiba','sechiba_init','','')
2758	END IF
2759
2760	!! 2. Read the default value that will be put into variable which are not in the restart file
2761	CALL ioget_expval(val_exp)
2762
2763	IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_init done '
2764
2765	END SUBROUTINE sechiba_init
2766
2767
2768	!! ==============================================================================================================================\n
2769	!! SUBROUTINE : sechiba_clear
2770	!!
2771	!>\BRIEF Deallocate memory of sechiba's variables
2772	!!
2773	!! DESCRIPTION : None
2774	!!
2775	!! RECENT CHANGE(S): None
2776	!!
2777	!! MAIN OUTPUT VARIABLE(S): None
2778	!!
2779	!! REFERENCE(S) : None
2780	!!
2781	!! FLOWCHART : None
2782	!! \n
2783	!_ ================================================================================================================================
2784
2785	SUBROUTINE sechiba_clear()
2786
2787	!! 1. Initialize first run
2788
2789	l_first_sechiba=.TRUE.
2790
2791	!! 2. Deallocate dynamic variables of sechiba
2792
2793	IF ( ALLOCATED (indexveg)) DEALLOCATE (indexveg)
2794	IF ( ALLOCATED (indexlai)) DEALLOCATE (indexlai)
2795	IF ( ALLOCATED (indexsoil)) DEALLOCATE (indexsoil)
2796	IF ( ALLOCATED (indexnobio)) DEALLOCATE (indexnobio)
2797	IF ( ALLOCATED (indexsnow)) DEALLOCATE (indexsnow)
2798	IF ( ALLOCATED (indexgrnd)) DEALLOCATE (indexgrnd)
2799	IF ( ALLOCATED (indexlayer)) DEALLOCATE (indexlayer)
2800	IF ( ALLOCATED (indexnslm)) DEALLOCATE (indexnslm)
2801	IF ( ALLOCATED (indexalb)) DEALLOCATE (indexalb)
2802	IF ( ALLOCATED (flood_res)) DEALLOCATE (flood_res)
2803	IF ( ALLOCATED (flood_frac)) DEALLOCATE (flood_frac)
2804	IF ( ALLOCATED (snow)) DEALLOCATE (snow)
2805	IF ( ALLOCATED (snow_age)) DEALLOCATE (snow_age)
2806	IF ( ALLOCATED (drysoil_frac)) DEALLOCATE (drysoil_frac)
2807	IF ( ALLOCATED (evap_bare_lim)) DEALLOCATE (evap_bare_lim)
2808	IF ( ALLOCATED (evap_bare_lim_ns)) DEALLOCATE (evap_bare_lim_ns)
2809	IF ( ALLOCATED (evapot)) DEALLOCATE (evapot)
2810	IF ( ALLOCATED (netrad)) DEALLOCATE (netrad)
2811	IF ( ALLOCATED (lwabs)) DEALLOCATE (lwabs)
2812	IF ( ALLOCATED (lwnet)) DEALLOCATE (lwnet)
2813	IF ( ALLOCATED (fluxsubli)) DEALLOCATE (fluxsubli)
2814	IF ( ALLOCATED (evapot_corr)) DEALLOCATE (evapot_corr)
2815	IF ( ALLOCATED (humrel)) DEALLOCATE (humrel)
2816	IF ( ALLOCATED (vegstress)) DEALLOCATE (vegstress)
2817	IF ( ALLOCATED (soiltile)) DEALLOCATE (soiltile)
2818	IF ( ALLOCATED (fraclut)) DEALLOCATE (fraclut)
2819	IF ( ALLOCATED (nwdFraclut)) DEALLOCATE (nwdFraclut)
2820	IF ( ALLOCATED (njsc)) DEALLOCATE (njsc)
2821	IF ( ALLOCATED (reinf_slope)) DEALLOCATE (reinf_slope)
2822	IF ( ALLOCATED (ks)) DEALLOCATE (ks)
2823	IF ( ALLOCATED (nvan)) DEALLOCATE (nvan)
2824	IF ( ALLOCATED (avan)) DEALLOCATE (avan)
2825	IF ( ALLOCATED (mcr)) DEALLOCATE (mcr)
2826	IF ( ALLOCATED (mcs)) DEALLOCATE (mcs)
2827	IF ( ALLOCATED (mcfc)) DEALLOCATE (mcfc)
2828	IF ( ALLOCATED (mcw)) DEALLOCATE (mcw)
2829	IF ( ALLOCATED (us)) DEALLOCATE (us)
2830	IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1)
2831	IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4)
2832	IF ( ALLOCATED (vbeta5)) DEALLOCATE (vbeta5)
2833	IF ( ALLOCATED (soilcap)) DEALLOCATE (soilcap)
2834	IF ( ALLOCATED (soilflx)) DEALLOCATE (soilflx)
2835	IF ( ALLOCATED (temp_sol)) DEALLOCATE (temp_sol)
2836	IF ( ALLOCATED (qsurf)) DEALLOCATE (qsurf)
2837	IF ( ALLOCATED (qsintveg)) DEALLOCATE (qsintveg)
2838	IF ( ALLOCATED (vbeta2)) DEALLOCATE (vbeta2)
2839	IF ( ALLOCATED (vbeta3)) DEALLOCATE (vbeta3)
2840	IF ( ALLOCATED (vbeta3pot)) DEALLOCATE (vbeta3pot)
2841	IF ( ALLOCATED (gsmean)) DEALLOCATE (gsmean)
2842	IF ( ALLOCATED (cimean)) DEALLOCATE (cimean)
2843	IF ( ALLOCATED (gpp)) DEALLOCATE (gpp)
2844	IF ( ALLOCATED (temp_growth)) DEALLOCATE (temp_growth)
2845	IF ( ALLOCATED (veget)) DEALLOCATE (veget)
2846	IF ( ALLOCATED (veget_max)) DEALLOCATE (veget_max)
2847	IF ( ALLOCATED (tot_bare_soil)) DEALLOCATE (tot_bare_soil)
2848	IF ( ALLOCATED (lai)) DEALLOCATE (lai)
2849	IF ( ALLOCATED (frac_age)) DEALLOCATE (frac_age)
2850	IF ( ALLOCATED (height)) DEALLOCATE (height)
2851	IF ( ALLOCATED (height_dom)) DEALLOCATE (height_dom)
2852	IF ( ALLOCATED (roughheight)) DEALLOCATE (roughheight)
2853	IF ( ALLOCATED (frac_nobio)) DEALLOCATE (frac_nobio)
2854	IF ( ALLOCATED (snow_nobio)) DEALLOCATE (snow_nobio)
2855	IF ( ALLOCATED (snow_nobio_age)) DEALLOCATE (snow_nobio_age)
2856	IF ( ALLOCATED (assim_param)) DEALLOCATE (assim_param)
2857	IF ( ALLOCATED (vevapflo)) DEALLOCATE (vevapflo)
2858	IF ( ALLOCATED (vevapsno)) DEALLOCATE (vevapsno)
2859	IF ( ALLOCATED (vevapnu)) DEALLOCATE (vevapnu)
2860	IF ( ALLOCATED (totfrac_nobio)) DEALLOCATE (totfrac_nobio)
2861	IF ( ALLOCATED (floodout)) DEALLOCATE (floodout)
2862	IF ( ALLOCATED (runoff)) DEALLOCATE (runoff)
2863	IF ( ALLOCATED (drainage)) DEALLOCATE (drainage)
2864	IF ( ALLOCATED (reinfiltration)) DEALLOCATE (reinfiltration)
2865	IF ( ALLOCATED (irrigation)) DEALLOCATE (irrigation)
2866	IF ( ALLOCATED (tot_melt)) DEALLOCATE (tot_melt)
2867	IF ( ALLOCATED (vbeta)) DEALLOCATE (vbeta)
2868	IF ( ALLOCATED (rau)) DEALLOCATE (rau)
2869	IF ( ALLOCATED (deadleaf_cover)) DEALLOCATE (deadleaf_cover)
2870	IF ( ALLOCATED (stempdiag)) DEALLOCATE (stempdiag)
2871	IF ( ALLOCATED (co2_to_bm)) DEALLOCATE (co2_to_bm)
2872	IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
2873	IF ( ALLOCATED (shumdiag_perma)) DEALLOCATE (shumdiag_perma)
2874	IF ( ALLOCATED (litterhumdiag)) DEALLOCATE (litterhumdiag)
2875	IF ( ALLOCATED (ptnlev1)) DEALLOCATE (ptnlev1)
2876	IF ( ALLOCATED (k_litt)) DEALLOCATE (k_litt)
2877	IF ( ALLOCATED (vevapwet)) DEALLOCATE (vevapwet)
2878	IF ( ALLOCATED (transpir)) DEALLOCATE (transpir)
2879	IF ( ALLOCATED (stressed)) DEALLOCATE (stressed)
2880	IF ( ALLOCATED (unstressed)) DEALLOCATE (unstressed)
2881	IF ( ALLOCATED (transpir_mod)) DEALLOCATE (transpir_mod)
2882	IF ( ALLOCATED (transpir_supply)) DEALLOCATE (transpir_supply)
2883	IF ( ALLOCATED (vir_transpir_supply)) DEALLOCATE (vir_transpir_supply)
2884	IF ( ALLOCATED (e_frac)) DEALLOCATE (e_frac)
2885	IF ( ALLOCATED (transpot)) DEALLOCATE (transpot)
2886	IF ( ALLOCATED (qsintmax)) DEALLOCATE (qsintmax)
2887	IF ( ALLOCATED (rveget)) DEALLOCATE (rveget)
2888	IF ( ALLOCATED (rstruct)) DEALLOCATE (rstruct)
2889	IF ( ALLOCATED (frac_snow_veg)) DEALLOCATE (frac_snow_veg)
2890	IF ( ALLOCATED (frac_snow_nobio)) DEALLOCATE (frac_snow_nobio)
2891	IF ( ALLOCATED (snowrho)) DEALLOCATE (snowrho)
2892	IF ( ALLOCATED (snowgrain)) DEALLOCATE (snowgrain)
2893	IF ( ALLOCATED (snowtemp)) DEALLOCATE (snowtemp)
2894	IF ( ALLOCATED (snowdz)) DEALLOCATE (snowdz)
2895	IF ( ALLOCATED (snowheat)) DEALLOCATE (snowheat)
2896	IF ( ALLOCATED (cgrnd_snow)) DEALLOCATE (cgrnd_snow)
2897	IF ( ALLOCATED (dgrnd_snow)) DEALLOCATE (dgrnd_snow)
2898	IF ( ALLOCATED (lambda_snow)) DEALLOCATE(lambda_snow)
2899	IF ( ALLOCATED (temp_sol_add)) DEALLOCATE(temp_sol_add)
2900	IF ( ALLOCATED (qsol_sat_new)) DEALLOCATE(qsol_sat_new)
2901	IF ( ALLOCATED (qair_new)) DEALLOCATE(qair_new)
2902	IF ( ALLOCATED (gtemp)) DEALLOCATE (gtemp)
2903	IF ( ALLOCATED (pgflux)) DEALLOCATE (pgflux)
2904	IF ( ALLOCATED (u_speed)) DEALLOCATE(u_speed)
2905	IF ( ALLOCATED (warnings)) DEALLOCATE (warnings)
2906	IF ( ALLOCATED (lai_per_level)) DEALLOCATE(lai_per_level)
2907	IF ( ALLOCATED (laieff_fit)) DEALLOCATE(laieff_fit)
2908	IF ( ALLOCATED (max_height_store)) DEALLOCATE(max_height_store)
2909	IF ( ALLOCATED (z_array_out)) DEALLOCATE(z_array_out)
2910	IF ( ALLOCATED (Light_Abs_Tot)) DEALLOCATE(Light_Abs_Tot)
2911	IF ( ALLOCATED (Light_Tran_Tot)) DEALLOCATE(Light_Tran_Tot)
2912	IF ( ALLOCATED (laieff_isotrop)) DEALLOCATE(laieff_isotrop)
2913	IF ( ALLOCATED (circ_class_biomass)) DEALLOCATE(circ_class_biomass)
2914	IF ( ALLOCATED (circ_class_n)) DEALLOCATE(circ_class_n)
2915	IF ( ALLOCATED (loss_gain)) DEALLOCATE(loss_gain)
2916	IF ( ALLOCATED (veget_max_new)) DEALLOCATE(veget_max_new)
2917	IF ( ALLOCATED (frac_nobio_new)) DEALLOCATE(frac_nobio_new)
2918	IF ( ALLOCATED (mc_layh)) DEALLOCATE (mc_layh)
2919	IF ( ALLOCATED (mcl_layh)) DEALLOCATE (mcl_layh)
2920	IF ( ALLOCATED (soilmoist)) DEALLOCATE (soilmoist)
2921	IF ( ALLOCATED (profile_vbeta3)) DEALLOCATE(profile_vbeta3)
2922	IF ( ALLOCATED (profile_rveget)) DEALLOCATE(profile_rveget)
2923	IF ( ALLOCATED (delta_c13_assim)) DEALLOCATE (delta_c13_assim)
2924	IF ( ALLOCATED (leaf_ci_out)) DEALLOCATE (leaf_ci_out)
2925	IF ( ALLOCATED (mcl_layh_s)) DEALLOCATE (mc_layh_s)
2926	IF ( ALLOCATED (mc_layh_s)) DEALLOCATE (mc_layh_s)
2927	IF ( ALLOCATED (soilmoist_s)) DEALLOCATE (soilmoist_s)
2928	IF ( ALLOCATED (tdeep)) DEALLOCATE (tdeep)
2929	IF ( ALLOCATED (hsdeep)) DEALLOCATE (hsdeep)
2930	IF ( ALLOCATED (heat_Zimov)) DEALLOCATE (heat_Zimov)
2931	IF ( ALLOCATED (sfluxCH4_deep)) DEALLOCATE (sfluxCH4_deep)
2932	IF ( ALLOCATED (sfluxCO2_deep)) DEALLOCATE (sfluxCO2_deep)
2933	IF ( ALLOCATED (som_total)) DEALLOCATE (som_total)
2934	IF ( ALLOCATED (altmax)) DEALLOCATE (altmax)
2935	IF ( ALLOCATED (swc)) DEALLOCATE(swc)
2936	IF ( ALLOCATED (ksoil)) DEALLOCATE(ksoil)
2937	IF ( ALLOCATED (root_profile)) DEALLOCATE (root_profile)
2938	IF ( ALLOCATED (coszang_noon)) DEALLOCATE (coszang_noon)
2939
2940	!! 3. Clear all allocated memory
2941
2942	CALL pft_parameters_clear
2943	CALL slowproc_clear
2944	CALL diffuco_clear
2945	CALL enerbil_clear
2946	CALL hydrol_clear
2947	CALL thermosoil_clear
2948	CALL condveg_clear
2949	CALL routing_clear
2950
2951	END SUBROUTINE sechiba_clear
2952
2953
2954	!! ==============================================================================================================================\n
2955	!! SUBROUTINE : sechiba_var_init
2956	!!
2957	!>\BRIEF Calculate air density as a function of air temperature and
2958	!! pressure for each terrestrial pixel.
2959	!!
2960	!! RECENT CHANGE(S): None
2961	!!
2962	!! MAIN OUTPUT VARIABLE(S): air density (::rau, kg m^{-3}).
2963	!!
2964	!! REFERENCE(S) : None
2965	!!
2966	!! FLOWCHART : None
2967	!! \n
2968	!_ ================================================================================================================================
2969
2970	SUBROUTINE sechiba_var_init (kjpindex, rau, pb, temp_air)
2971
2972	!! 0.1 Input variables
2973
2974	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
2975	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
2976	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
2977
2978	!! 0.2 Output variables
2979
2980	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rau !! Air density @tex $(kg m^{-3})$ @endtex
2981
2982	!! 0.3 Modified variables
2983
2984	!! 0.4 Local variables
2985
2986	INTEGER(i_std) :: ji !! Indices (unitless)
2987	!_ ================================================================================================================================
2988
2989	!! 1. Calculate intial air density (::rau)
2990
2991	DO ji = 1,kjpindex
2992	rau(ji) = pa_par_hpa * pb(ji) / (cte_molr*temp_air(ji))
2993	END DO
2994
2995	IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_var_init done '
2996
2997	END SUBROUTINE sechiba_var_init
2998
2999
3000	!! ==============================================================================================================================\n
3001	!! SUBROUTINE : sechiba_end
3002	!!
3003	!>\BRIEF Swap old for newly calculated soil temperature.
3004	!!
3005	!! RECENT CHANGE(S): None
3006	!!
3007	!! MAIN OUTPUT VARIABLE(S): soil temperature (::temp_sol; K)
3008	!!
3009	!! REFERENCE(S) : None
3010	!!
3011	!! FLOWCHART : None
3012	!! \n
3013	!! ================================================================================================================================
3014
3015	SUBROUTINE sechiba_end (kjpindex, temp_sol_new, temp_sol)
3016
3017
3018	!! 0.1 Input variables
3019
3020	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
3021	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature (K)
3022
3023	!! 0.2 Output variables
3024	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol !! Soil temperature (K)
3025
3026	!_ ================================================================================================================================
3027
3028	!! 1. Swap temperature
3029
3030	temp_sol(:) = temp_sol_new(:)
3031
3032	IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_end done '
3033
3034	END SUBROUTINE sechiba_end
3035
3036
3037	!! ==============================================================================================================================\n
3038	!! SUBROUTINE : sechiba_interface_orchidee_inca
3039	!!
3040	!>\BRIEF make the interface between surface and atmospheric chemistry
3041	!!
3042	!! DESCRIPTION : This subroutine is called from INCA, the atmospheric chemistry model. It is used to transfer variables from ORCHIDEE to INCA.
3043	!!
3044	!! RECENT CHANGE(S): move from chemistry module to be more generic (feb - 2017)
3045	!!
3046	!! MAIN OUTPUT VARIABLE(S): emission COV to be transport by orchidee to inca in fields_out array
3047	!!
3048	!! REFERENCE(S) : None
3049	!!
3050	!! FLOWCHART : None
3051	!! \n
3052	!! ================================================================================================================================
3053	SUBROUTINE sechiba_interface_orchidee_inca( &
3054	nvm_out, veget_max_out, veget_frac_out, lai_out, snow_out, &
3055	field_out_names, fields_out, field_in_names, fields_in)
3056
3057
3058	INTEGER, INTENT(out) :: nvm_out !! Number of vegetation types
3059	REAL(r_std), DIMENSION (:,:), INTENT(out) :: veget_max_out !! Max. fraction of vegetation type (LAI -> infty)
3060	REAL(r_std), DIMENSION (:,:), INTENT(out) :: veget_frac_out !! Fraction of vegetation type (unitless, 0-1)
3061	REAL(r_std), DIMENSION (:,:), INTENT(out) :: lai_out !! Surface foliere
3062	REAL(r_std), DIMENSION (:) , INTENT(out) :: snow_out !! Snow mass [Kg/m^2]
3063
3064	!
3065	! Optional arguments
3066	!
3067	! Names and fields for emission variables : to be transport by Orchidee to Inca
3068	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_out_names
3069	REAL(r_std),DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: fields_out
3070	!
3071	! Names and fields for deposit variables : to be transport from chemistry model by INCA to ORCHIDEE.
3072	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_in_names
3073	REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(IN) :: fields_in
3074
3075
3076	! Variables always transmitted from sechiba to inca
3077	nvm_out = nvm
3078	veget_max_out(:,:) = veget_max(:,:)
3079	veget_frac_out(:,:) = veget(:,:)
3080	lai_out(:,:) = lai(:,:)
3081	snow_out(:) = snow(:)
3082
3083	! Call chemistry_flux_interface if at least one of variables field_out_names or
3084	! field_in_names is present in the argument list of sechiba_interface_orchidee_inca when called from inca.
3085	IF (PRESENT(field_out_names) .AND. .NOT. PRESENT(field_in_names)) THEN
3086	CALL chemistry_flux_interface(field_out_names=field_out_names, fields_out=fields_out)
3087	ELSE IF (.NOT. PRESENT(field_out_names) .AND. PRESENT(field_in_names)) THEN
3088	CALL chemistry_flux_interface(field_in_names=field_in_names, fields_in=fields_in)
3089	ELSE IF (PRESENT(field_out_names) .AND. PRESENT(field_in_names)) THEN
3090	CALL chemistry_flux_interface(field_out_names=field_out_names, fields_out=fields_out, &
3091	field_in_names=field_in_names, fields_in=fields_in)
3092	ENDIF
3093
3094	END SUBROUTINE sechiba_interface_orchidee_inca
3095
3096	END MODULE sechiba
3097

Note: See TracBrowser for help on using the repository browser.

Download in other formats: