Context Navigation

source: branches/publications/ORCHIDEE_CN_CAN_r5698/src_sechiba/sechiba.f90 @ 7346

Last change on this file since 7346 was 5691, checked in by sebastiaan.luyssaert, 6 years ago
DEV: tested with 13, 37 and 64 PFTs with LCC on different pixels. Some configuration run for 20 years on a given pixel, other crash on another pixel. There is a mass balance problem in sapiens_lcc (ticket #482). This commit fixes a problem with PFT1 in littercalc. This PFT is now fully integrated in LCC and subsequent litter and soil dynamics. veget_max was changed in veget_cov_max where appropriate, a typo in enerbil was corrected.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 174.8 KB

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

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

Download in other formats: