Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90 @ 8221

Last change on this file since 8221 was 8221, checked in by anne.cozic, 8 months ago
update sechiba_interface_orchidee_inca to allow compatibility between all versions of Orchidee and a single one for Inca
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 117.7 KB

Line
1	! ==============================================================================================================================\n
2	! MODULE : sechiba
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Structures the calculation of atmospheric and hydrological
10	!! variables by calling diffuco_main, enerbil_main, hydrol_main,
11	!! condveg_main and thermosoil_main. Note that sechiba_main
12	!! calls slowproc_main and thus indirectly calculates the biogeochemical
13	!! processes as well.
14	!!
15	!!\n DESCRIPTION : :: shumdiag, :: litterhumdiag and :: stempdiag :: ftempdiag are not
16	!! saved in the restart file because at the first time step because they
17	!! are recalculated. However, they must be saved as they are in slowproc
18	!! which is called before the modules which calculate them.
19	!!
20	!! RECENT CHANGE(S): November 2020: It is possible to define soil hydraulic parameters from maps,
21	!! as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes).
22	!! Here, it leads to declare and allocate global variables.
23	!!
24	!! REFERENCE(S) : None
25	!!
26	!! SVN :
27	!! $HeadURL$
28	!! $Date$
29	!! $Revision$
30	!! \n
31	!_ ================================================================================================================================
32
33	MODULE sechiba
34
35	USE ioipsl
36	USE xios_orchidee
37	USE constantes
38	USE time, ONLY : one_day, dt_sechiba
39	USE constantes_soil
40	USE pft_parameters
41	USE grid
42	USE diffuco
43	USE condveg
44	USE enerbil
45	USE hydrol
46	USE thermosoil
47	USE sechiba_io_p
48	USE slowproc
49	USE routing_wrapper
50	USE ioipsl_para
51	USE chemistry
52
53	IMPLICIT NONE
54
55	PRIVATE
56	PUBLIC sechiba_main, sechiba_initialize, sechiba_clear, sechiba_interface_orchidee_inca, sechiba_xios_initialize
57
58	INTEGER(i_std), SAVE :: printlev_loc !! local printlev for this module
59	!$OMP THREADPRIVATE(printlev_loc)
60	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexveg !! indexing array for the 3D fields of vegetation
61	!$OMP THREADPRIVATE(indexveg)
62	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlai !! indexing array for the 3D fields of vegetation
63	!$OMP THREADPRIVATE(indexlai)
64	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnobio !! indexing array for the 3D fields of other surfaces (ice,
65	!! lakes, ...)
66	!$OMP THREADPRIVATE(indexnobio)
67	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsoil !! indexing array for the 3D fields of soil types (kjpindex*nstm)
68	!$OMP THREADPRIVATE(indexsoil)
69	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexgrnd !! indexing array for the 3D ground heat profiles (kjpindex*ngrnd)
70	!$OMP THREADPRIVATE(indexgrnd)
71	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlayer !! indexing array for the 3D fields of soil layers in CWRR (kjpindex*nslm)
72	!$OMP THREADPRIVATE(indexlayer)
73	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnslm !! indexing array for the 3D fields of diagnostic soil layers (kjpindex*nslm)
74	!$OMP THREADPRIVATE(indexnslm)
75	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexalb !! indexing array for the 2 fields of albedo
76	!$OMP THREADPRIVATE(indexalb)
77	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsnow !! indexing array for the 3D fields snow layers
78	!$OMP THREADPRIVATE(indexsnow)
79	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget !! Fraction of vegetation type (unitless, 0-1)
80	!$OMP THREADPRIVATE(veget)
81	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget_max !! Max. fraction of vegetation type (LAI -> infty, unitless)
82	!$OMP THREADPRIVATE(veget_max)
83	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_bare_soil !! Total evaporating bare soil fraction
84	!$OMP THREADPRIVATE(tot_bare_soil)
85	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: height !! Vegetation Height (m)
86	!$OMP THREADPRIVATE(height)
87	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
88	!! (unitless, 0-1)
89	!$OMP THREADPRIVATE(totfrac_nobio)
90	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: floodout !! Flow out of floodplains from hydrol
91	!$OMP THREADPRIVATE(floodout)
92	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: runoff !! Surface runoff calculated by hydrol
93	!! @tex $(kg m^{-2})$ @endtex
94	!$OMP THREADPRIVATE(runoff)
95	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: drainage !! Deep drainage calculatedd by hydrol
96	!! @tex $(kg m^{-2})$ @endtex
97	!$OMP THREADPRIVATE(drainage)
98	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: returnflow !! Water flow from lakes and swamps which returns to
99	!! the grid box @tex $(kg m^{-2})$ @endtex
100	!$OMP THREADPRIVATE(returnflow)
101	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: reinfiltration !! Routed water which returns into the soil
102	!$OMP THREADPRIVATE(reinfiltration)
103	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: irrigation !! Irrigation flux taken from the routing reservoirs and
104	!! being put into the upper layers of the soil
105	!! @tex $(kg m^{-2})$ @endtex
106	!$OMP THREADPRIVATE(irrigation)
107	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: emis !! Surface emissivity (unitless)
108	!$OMP THREADPRIVATE(emis)
109	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z0h !! Surface roughness for heat (m)
110	!$OMP THREADPRIVATE(z0h)
111	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z0m !! Surface roughness for momentum (m)
112	!$OMP THREADPRIVATE(z0m)
113	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: roughheight !! Effective height for roughness (m)
114	!$OMP THREADPRIVATE(roughheight)
115	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: reinf_slope !! slope coefficient (reinfiltration)
116	!$OMP THREADPRIVATE(reinf_slope)
117	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: reinf_slope_soil !! slope coefficient (reinfiltration) per soil tile
118	!$OMP THREADPRIVATE(reinf_slope_soil)
119
120
121	!salma: adding soil hydraulic params
122	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: ks !! Saturated soil conductivity (mm d^{-1})
123	!$OMP THREADPRIVATE(ks)
124	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: nvan !! Van Genushten n parameter (unitless)
125	!$OMP THREADPRIVATE(nvan)
126	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: avan !! Van Genushten alpha parameter (mm ^{-1})
127	!$OMP THREADPRIVATE(avan)
128	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcr !! Residual soil moisture (m^{3} m^{-3})
129	!$OMP THREADPRIVATE(mcr)
130	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcs !! Saturated soil moisture (m^{3} m^{-3})
131	!$OMP THREADPRIVATE(mcs)
132	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcfc !! Volumetric water content at field capacity (m^{3} m^{-3})
133	!$OMP THREADPRIVATE(mcfc)
134	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcw !! Volumetric water content at wilting point (m^{3} m^{-3})
135	!$OMP THREADPRIVATE(mcw)
136	!end salma:adding soil hydraulic params
137
138
139	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shumdiag !! Mean relative soil moisture in the different levels used
140	!! by thermosoil.f90 (unitless, 0-1)
141	!$OMP THREADPRIVATE(shumdiag)
142	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shumdiag_perma !! Saturation degree of the soil
143	!$OMP THREADPRIVATE(shumdiag_perma)
144	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: k_litt !! litter cond.
145	!$OMP THREADPRIVATE(k_litt)
146	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: litterhumdiag !! Litter dryness factor (unitless, 0-1)
147	!$OMP THREADPRIVATE(litterhumdiag)
148	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: stempdiag !! Temperature which controls canopy evolution (K)
149	!$OMP THREADPRIVATE(stempdiag)
150	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ftempdiag !! Temperature over the full soil column for river temperature (K)
151	!$OMP THREADPRIVATE(ftempdiag)
152	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: qsintveg !! Water on vegetation due to interception
153	!! @tex $(kg m^{-2})$ @endtex
154	!$OMP THREADPRIVATE(qsintveg)
155	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vbeta2 !! Interception resistance (unitless,0-1)
156	!$OMP THREADPRIVATE(vbeta2)
157	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vbeta3 !! Vegetation resistance (unitless,0-1)
158	!$OMP THREADPRIVATE(vbeta3)
159	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vbeta3pot !! Potential vegetation resistance
160	!$OMP THREADPRIVATE(vbeta3pot)
161	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gsmean !! Mean stomatal conductance for CO2 (mol m-2 s-1)
162	!$OMP THREADPRIVATE(gsmean)
163	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cimean !! STOMATE: mean intercellular CO2 concentration (ppm)
164	!$OMP THREADPRIVATE(cimean)
165	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vevapwet !! Interception loss over each PFT
166	!! @tex $(kg m^{-2} days^{-1})$ @endtex
167	!$OMP THREADPRIVATE(vevapwet)
168	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: transpir !! Transpiration @tex $(kg m^{-2} days^{-1})$ @endtex
169	!$OMP THREADPRIVATE(transpir)
170	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: transpot !! Potential Transpiration (needed for irrigation)
171	!$OMP THREADPRIVATE(transpot)
172	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: qsintmax !! Maximum amount of water in the canopy interception
173	!! reservoir @tex $(kg m^{-2})$ @endtex
174	!$OMP THREADPRIVATE(qsintmax)
175	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: rveget !! Surface resistance for the vegetation
176	!! @tex $(s m^{-1})$ @endtex
177	!$OMP THREADPRIVATE(rveget)
178	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: rstruct !! Vegetation structural resistance
179	!$OMP THREADPRIVATE(rstruct)
180	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: snow_nobio !! Snow mass of non-vegetative surfaces
181	!! @tex $(kg m^{-2})$ @endtex
182	!$OMP THREADPRIVATE(snow_nobio)
183	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: snow_nobio_age !! Snow age on non-vegetative surfaces (days)
184	!$OMP THREADPRIVATE(snow_nobio_age)
185	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: frac_nobio !! Fraction of non-vegetative surfaces (continental ice,
186	!! lakes, ...) (unitless, 0-1)
187	!$OMP THREADPRIVATE(frac_nobio)
188	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: assim_param !! min+max+opt temps, vcmax, vjmax for photosynthesis
189	!$OMP THREADPRIVATE(assim_param)
190	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: lai !! Surface foliaire
191	!$OMP THREADPRIVATE(lai)
192	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gpp !! STOMATE: GPP. gC/m**2 of total area
193	!$OMP THREADPRIVATE(gpp)
194	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: temp_growth !! Growth temperature (C) - Is equal to t2m_month
195	!$OMP THREADPRIVATE(temp_growth)
196	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: humrel !! Relative humidity
197	!$OMP THREADPRIVATE(humrel)
198	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vegstress !! Vegetation moisture stress (only for vegetation growth)
199	!$OMP THREADPRIVATE(vegstress)
200	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: frac_age !! Age efficacity from STOMATE for isoprene
201	!$OMP THREADPRIVATE(frac_age)
202	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: soiltile !! Fraction of each soil tile (0-1, unitless)
203	!$OMP THREADPRIVATE(soiltile)
204	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: fraclut !! Fraction of each landuse tile (0-1, unitless)
205	!$OMP THREADPRIVATE(fraclut)
206	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: nwdFraclut !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless)
207	!$OMP THREADPRIVATE(nwdFraclut)
208	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
209	!$OMP THREADPRIVATE(njsc)
210	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta1 !! Snow resistance
211	!$OMP THREADPRIVATE(vbeta1)
212	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta4 !! Bare soil resistance
213	!$OMP THREADPRIVATE(vbeta4)
214	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta5 !! Floodplains resistance
215	!$OMP THREADPRIVATE(vbeta5)
216	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: soilcap !!
217	!$OMP THREADPRIVATE(soilcap)
218	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: soilflx !!
219	!$OMP THREADPRIVATE(soilflx)
220	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol !! Soil temperature
221	!$OMP THREADPRIVATE(temp_sol)
222	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsurf !! near soil air moisture
223	!$OMP THREADPRIVATE(qsurf)
224	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: flood_res !! flood reservoir estimate
225	!$OMP THREADPRIVATE(flood_res)
226	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: flood_frac !! flooded fraction
227	!$OMP THREADPRIVATE(flood_frac)
228	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow !! Snow mass [Kg/m^2]
229	!$OMP THREADPRIVATE(snow)
230	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_age !! Snow age @tex ($d$) @endtex
231	!$OMP THREADPRIVATE(snow_age)
232	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: drysoil_frac !! Fraction of visibly (albedo) Dry soil (Between 0 and 1)
233	!$OMP THREADPRIVATE(drysoil_frac)
234	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: evap_bare_lim !! Bare soil stress
235	!$OMP THREADPRIVATE(evap_bare_lim)
236	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: evap_bare_lim_ns !! Bare soil stress
237	!$OMP THREADPRIVATE(evap_bare_lim_ns)
238	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: co2_flux !! CO2 flux (gC/m**2 of average ground/one_day)
239	!$OMP THREADPRIVATE(co2_flux)
240	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: evapot !! Soil Potential Evaporation
241	!$OMP THREADPRIVATE(evapot)
242	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: evapot_corr !! Soil Potential Evaporation Correction (Milly 1992)
243	!$OMP THREADPRIVATE(evapot_corr)
244	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vevapflo !! Floodplains evaporation
245	!$OMP THREADPRIVATE(vevapflo)
246	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vevapsno !! Snow evaporation
247	!$OMP THREADPRIVATE(vevapsno)
248	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vevapnu !! Bare soil evaporation
249	!$OMP THREADPRIVATE(vevapnu)
250	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_melt !! Total melt
251	!$OMP THREADPRIVATE(tot_melt)
252	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta !! Resistance coefficient
253	!$OMP THREADPRIVATE(vbeta)
254	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: rau !! Density
255	!$OMP THREADPRIVATE(rau)
256	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: deadleaf_cover !! Fraction of soil covered by dead leaves
257	!$OMP THREADPRIVATE(deadleaf_cover)
258	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: ptnlev1 !! 1st level Different levels soil temperature
259	!$OMP THREADPRIVATE(ptnlev1)
260	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mc_layh !! Volumetric soil moisture for each layer in hydrol(liquid + ice) (m3/m3)
261	!$OMP THREADPRIVATE(mc_layh)
262	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mcl_layh !! Volumetric soil moisture for each layer in hydrol(liquid) (m3/m3)
263	!$OMP THREADPRIVATE(mcl_layh)
264	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: soilmoist !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
265	!$OMP THREADPRIVATE(soilmoist)
266
267	LOGICAL, SAVE :: l_first_sechiba = .TRUE. !! Flag controlling the intialisation (true/false)
268	!$OMP THREADPRIVATE(l_first_sechiba)
269
270	! Variables related to snow processes calculations
271
272	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: frac_snow_veg !! Snow cover fraction on vegetation (unitless)
273	!$OMP THREADPRIVATE(frac_snow_veg)
274	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: frac_snow_nobio !! Snow cover fraction on continental ice, lakes, etc (unitless)
275	!$OMP THREADPRIVATE(frac_snow_nobio)
276	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: snowrho !! snow density for each layer
277	!$OMP THREADPRIVATE(snowrho)
278	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: snowheat !! snow heat content for each layer (J/m2)
279	!$OMP THREADPRIVATE(snowheat)
280	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: snowgrain !! snow grain size (m)
281	!$OMP THREADPRIVATE(snowgrain)
282	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: snowtemp !! snow temperature profile (K)
283	!$OMP THREADPRIVATE(snowtemp)
284	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: snowdz !! snow layer thickness (m)
285	!$OMP THREADPRIVATE(snowdz)
286	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: gtemp !! soil surface temperature
287	!$OMP THREADPRIVATE(gtemp)
288	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: pgflux !! net energy into snow pack
289	!$OMP THREADPRIVATE(pgflux)
290	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
291	!$OMP THREADPRIVATE(cgrnd_snow)
292	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
293	!$OMP THREADPRIVATE(dgrnd_snow)
294	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
295	!! from the first and second snow layers
296	!$OMP THREADPRIVATE(lambda_snow)
297	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
298	!$OMP THREADPRIVATE(temp_sol_add)
299	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: root_deficit !! water deficit to reach IRRIGATION SOIL MOISTURE TARGET
300	!$OMP THREADPRIVATE(root_deficit)
301	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: irrig_frac !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
302	!$OMP THREADPRIVATE(irrig_frac)
303	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: irrigated_next !! Dynamic irrig. area, calculated in slowproc and passed to routing
304	!! @tex $(m^{-2})$ @endtex
305	!$OMP THREADPRIVATE(irrigated_next)
306	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: fraction_aeirrig_sw !! Fraction of area equipped for irrigation from surface water, of irrig_frac
307	!! 1.0 here corresponds to irrig_frac, not grid cell
308	!$OMP THREADPRIVATE(fraction_aeirrig_sw)
309
310	CONTAINS
311
312
313	!! =============================================================================================================================
314	!! SUBROUTINE: sechiba_xios_initialize
315	!!
316	!>\BRIEF Initialize xios dependant defintion before closing context defintion
317	!!
318	!! DESCRIPTION: Initialize xios dependant defintion before closing context defintion
319	!!
320	!! \n
321	!_ ==============================================================================================================================
322
323	SUBROUTINE sechiba_xios_initialize
324	LOGICAL :: lerr
325
326	IF (xios_orchidee_ok) THEN
327	lerr=xios_orchidee_setvar('min_sechiba',min_sechiba)
328	CALL slowproc_xios_initialize
329	CALL condveg_xios_initialize
330	CALL chemistry_xios_initialize
331	CALL thermosoil_xios_initialize
332	CALL routing_wrapper_xios_initialize
333	END IF
334	IF (printlev_loc>=3) WRITE(numout,*) 'End sechiba_xios_initialize'
335
336	END SUBROUTINE sechiba_xios_initialize
337
338
339
340
341	!! =============================================================================================================================
342	!! SUBROUTINE: sechiba_initialize
343	!!
344	!>\BRIEF Initialize all prinicipal modules by calling their "_initialize" subroutines
345	!!
346	!! DESCRIPTION: Initialize all prinicipal modules by calling their "_initialize" subroutines
347	!!
348	!! \n
349	!_ ==============================================================================================================================
350
351	SUBROUTINE sechiba_initialize( &
352	kjit, kjpij, kjpindex, index, &
353	lalo, contfrac, neighbours, resolution, zlev, &
354	u, v, qair, temp_air, &
355	petAcoef, peqAcoef, petBcoef, peqBcoef, &
356	precip_rain, precip_snow, lwdown, swnet, swdown, &
357	pb, rest_id, hist_id, hist2_id, &
358	rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
359	coastalflow, riverflow, tsol_rad, vevapp, qsurf_out, &
360	z0m_out, z0h_out, albedo, fluxsens, fluxlat, emis_out, &
361	temp_sol_new, tq_cdrag)
362
363	!! 0.1 Input variables
364	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
365	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
366	!! (unitless)
367	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
368	!! (unitless)
369	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
370	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (unitless)
371	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier (unitless)
372	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
373	!! (unitless)
374	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _History_ file identifier
375	!! (unitless)
376	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file
377	!! identifier (unitless)
378	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographic coordinates (latitude,longitude)
379	!! for grid cells (degrees)
380	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid
381	!! (unitless, 0-1)
382	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the pixels on the map.
383	!! Sechiba uses a reduced grid excluding oceans
384	!! ::index contains the indices of the
385	!! terrestrial pixels only! (unitless)
386	INTEGER(i_std), DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours !! Neighboring grid points if land!(unitless)
387	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size in x and y of the grid (m)
388
389	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed in direction u
390	!! @tex $(m.s^{-1})$ @endtex
391	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed in direction v
392	!! @tex $(m.s^{-1})$ @endtex
393	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
394	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
395	!! @tex $(kg kg^{-1})$ @endtex
396	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
397	!! @tex $(kg m^{-2})$ @endtex
398	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
399	!! @tex $(kg m^{-2})$ @endtex
400	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux
401	!! @tex $(W m^{-2})$ @endtex
402	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux
403	!! @tex $(W m^{-2})$ @endtex
404	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux
405	!! @tex $(W m^{-2})$ @endtex
406	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
407	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! Coefficients A for T from the Planetary
408	!! Boundary Layer
409	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! Coefficients A for q from the Planetary
410	!! Boundary Layer
411	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! Coefficients B for T from the Planetary
412	!! Boundary Layer
413	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! Coefficients B for q from the Planetary
414	!! Boundary Layer
415	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
416
417
418	!! 0.2 Output variables
419	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: coastalflow !! Outflow on coastal points by small basins.
420	!! This is the water which flows in a disperse
421	!! way into the ocean
422	!! @tex $(kg dt_routing^{-1})$ @endtex
423	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: riverflow !! Outflow of the major rivers.
424	!! The flux will be located on the continental
425	!! grid but this should be a coastal point
426	!! @tex $(kg dt_routing^{-1})$ @endtex
427	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature
428	!! @tex $(W m^{-2})$ @endtex
429	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation
430	!! @tex $(kg m^{-2} days^{-1})$ @endtex
431
432	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf_out !! Surface specific humidity
433	!! @tex $(kg kg^{-1})$ @endtex
434	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m_out !! Surface roughness momentum (output diagnostic, m)
435	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h_out !! Surface roughness heat (output diagnostic, m)
436	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Surface albedo for visible and near-infrared
437	!! (unitless, 0-1)
438	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux
439	!! @tex $(W m^{-2})$ @endtex
440	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux
441	!! @tex $(W m^{-2})$ @endtex
442	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis_out !! Emissivity (unitless)
443	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New ground temperature (K)
444
445	!! 0.3 Modified
446	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: tq_cdrag !! Surface drag coefficient (-)
447
448	!! 0.4 Local variables
449	INTEGER(i_std) :: ji, jv !! Index (unitless)
450	REAL(r_std), DIMENSION(kjpindex) :: zmaxh_glo !! 2D field of constant soil depth (zmaxh) (m)
451	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
452
453	!_ ================================================================================================================================
454
455	IF (printlev>=3) WRITE(numout,*) 'Start sechiba_initialize'
456
457	!! 1. Initialize variables on first call
458
459	!! 1.1 Initialize most of sechiba's variables
460	CALL sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)
461
462	!! 1.3 Initialize stomate's variables
463
464	CALL slowproc_initialize (kjit, kjpij, kjpindex, &
465	rest_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
466	index, indexveg, lalo, neighbours, &
467	resolution, contfrac, temp_air, &
468	soiltile, reinf_slope, ks, nvan, &
469	avan, mcr, mcs, mcfc, &
470	mcw, deadleaf_cover, assim_param, &
471	lai, frac_age, height, veget, &
472	frac_nobio, njsc, veget_max, fraclut, &
473	nwdfraclut, tot_bare_soil,totfrac_nobio, qsintmax, &
474	temp_growth, irrigated_next, irrig_frac, fraction_aeirrig_sw, &
475	reinf_slope_soil)
476
477	!! 1.4 Initialize diffusion coefficients
478	CALL diffuco_initialize (kjit, kjpindex, index, &
479	rest_id, lalo, neighbours, resolution, &
480	rstruct, tq_cdrag)
481
482	!! 1.5 Initialize remaining variables of energy budget
483	CALL enerbil_initialize (kjit, kjpindex, index, rest_id, &
484	qair, &
485	temp_sol, temp_sol_new, tsol_rad, &
486	evapot, evapot_corr, qsurf, fluxsens, &
487	fluxlat, vevapp )
488
489
490	!! 1.7 Initialize remaining hydrological variables
491	CALL hydrol_initialize (ks, nvan, avan, mcr, mcs, mcfc, mcw, &
492	kjit, kjpindex, index, rest_id, &
493	njsc, soiltile, veget, veget_max, &
494	humrel, vegstress, drysoil_frac, &
495	shumdiag_perma, qsintveg, &
496	evap_bare_lim, evap_bare_lim_ns, snow, snow_age, snow_nobio, &
497	snow_nobio_age, snowrho, snowtemp, snowgrain, &
498	snowdz, snowheat, &
499	mc_layh, mcl_layh, soilmoist)
500
501
502	!! 1.9 Initialize surface parameters (emissivity, albedo and roughness)
503	CALL condveg_initialize (kjit, kjpindex, index, rest_id, &
504	lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
505	zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
506	drysoil_frac, height, snowdz,snowrho, tot_bare_soil, &
507	temp_air, pb, u, v, lai, &
508	emis, albedo, z0m, z0h, roughheight, &
509	frac_snow_veg,frac_snow_nobio)
510
511
512	!! 1.10 Initialization of soil thermodynamics
513	CALL thermosoil_initialize (kjit, kjpindex, rest_id, mcs, &
514	temp_sol_new, snow, shumdiag_perma, &
515	soilcap, soilflx, stempdiag, ftempdiag, &
516	gtemp, &
517	mc_layh, mcl_layh, soilmoist, njsc , &
518	frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
519	snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb, &
520	ptnlev1)
521
522
523	!! 1.12 Initialize river routing
524	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
525
526	!! 1.12.1 Initialize river routing
527	CALL routing_wrapper_initialize( &
528	kjit, kjpindex, index, &
529	rest_id, hist_id, hist2_id, lalo, &
530	neighbours, resolution, contfrac, stempdiag, ftempdiag, &
531	soiltile, irrig_frac, veget_max, irrigated_next, &
532	returnflow, reinfiltration, irrigation, riverflow, &
533	coastalflow, flood_frac, flood_res )
534	ELSE
535	!! 1.12.2 No routing, set variables to zero
536	riverflow(:) = zero
537	coastalflow(:) = zero
538	returnflow(:) = zero
539	reinfiltration(:) = zero
540	irrigation(:) = zero
541	flood_frac(:) = zero
542	flood_res(:) = zero
543	ENDIF
544
545	!! 1.13 Write internal variables to output fields
546	z0m_out(:) = z0m(:)
547	z0h_out(:) = z0h(:)
548	emis_out(:) = emis(:)
549	qsurf_out(:) = qsurf(:)
550
551	!! 2. Output variables only once
552	zmaxh_glo(:) = zmaxh
553	CALL xios_orchidee_send_field("zmaxh",zmaxh_glo)
554
555	IF (printlev_loc>=3) WRITE(numout,*) 'sechiba_initialize done'
556
557	END SUBROUTINE sechiba_initialize
558
559	!! ==============================================================================================================================\n
560	!! SUBROUTINE : sechiba_main
561	!!
562	!>\BRIEF Main routine for the sechiba module performing three functions:
563	!! calculating temporal evolution of all variables and preparation of output and
564	!! restart files (during the last call only)
565	!!
566	!!\n DESCRIPTION : Main routine for the sechiba module.
567	!! One time step evolution consists of:
568	!! - call sechiba_var_init to do some initialization,
569	!! - call slowproc_main to do some daily calculations
570	!! - call diffuco_main for diffusion coefficient calculation,
571	!! - call enerbil_main for energy budget calculation,
572	!! - call hydrol_main for hydrologic processes calculation,
573	!! - call condveg_main for surface conditions such as roughness, albedo, and emmisivity,
574	!! - call thermosoil_main for soil thermodynamic calculation,
575	!! - call sechiba_end to swap previous to new fields.
576	!!
577	!! RECENT CHANGE(S): None
578	!!
579	!! MAIN OUTPUT VARIABLE(S): Hydrological variables (:: coastalflow and :: riverflow),
580	!! components of the energy budget (:: tsol_rad, :: vevapp, :: fluxsens,
581	!! :: temp_sol_new and :: fluxlat), surface characteristics (:: z0_out, :: emis_out,
582	!! :: tq_cdrag and :: albedo) and land use related CO2 fluxes (:: netco2flux and
583	!! :: fco2_lu, :: fco2_wh, ::fco2_ha)
584	!!
585	!! REFERENCE(S) :
586	!!
587	!! FLOWCHART :
588	!! \latexonly
589	!! \includegraphics[scale = 0.5]{sechibamainflow.png}
590	!! \endlatexonly
591	!! \n
592	!_ ================================================================================================================================
593
594	SUBROUTINE sechiba_main (kjit, kjpij, kjpindex, index, &
595	& ldrestart_read, ldrestart_write, &
596	& lalo, contfrac, neighbours, resolution,&
597	& zlev, u, v, qair, temp_air, epot_air, ccanopy, &
598	& tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
599	& precip_rain, precip_snow, lwdown, swnet, swdown, coszang, pb, &
600	& vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
601	& netco2flux, fco2_lu, fco2_wh, fco2_ha, &
602	& tsol_rad, temp_sol_new, qsurf_out, albedo, emis_out, z0m_out, z0h_out,&
603	& veget_out, lai_out, height_out, &
604	& rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
605
606	!! 0.1 Input variables
607
608	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
609	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
610	!! (unitless)
611	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
612	!! (unitless)
613	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
614	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (unitless)
615	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier (unitless)
616	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
617	!! (unitless)
618	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _History_ file identifier
619	!! (unitless)
620	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file
621	!! identifier (unitless)
622	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
623	!! (true/false)
624	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
625	!! (true/false)
626	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographic coordinates (latitude,longitude)
627	!! for grid cells (degrees)
628	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid
629	!! (unitless, 0-1)
630	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the pixels on the map.
631	!! Sechiba uses a reduced grid excluding oceans
632	!! ::index contains the indices of the
633	!! terrestrial pixels only! (unitless)
634	INTEGER(i_std), DIMENSION(kjpindex,NbNeighb), INTENT(in) :: neighbours !! Neighboring grid points if land!(unitless)
635	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size in x and y of the grid (m)
636
637	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed in direction u
638	!! @tex $(m.s^{-1})$ @endtex
639	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed in direction v
640	!! @tex $(m.s^{-1})$ @endtex
641	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
642	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
643	!! @tex $(kg kg^{-1})$ @endtex
644	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
645	!! @tex $(kg m^{-2})$ @endtex
646	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
647	!! @tex $(kg m^{-2})$ @endtex
648	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux
649	!! @tex $(W m^{-2})$ @endtex
650	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: coszang !! Cosine of the solar zenith angle (unitless)
651	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux
652	!! @tex $(W m^{-2})$ @endtex
653	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux
654	!! @tex $(W m^{-2})$ @endtex
655	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
656	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (??J)
657	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration in the canopy (ppm)
658	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! Coefficients A for T from the Planetary
659	!! Boundary Layer
660	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! Coefficients A for q from the Planetary
661	!! Boundary Layer
662	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! Coefficients B for T from the Planetary
663	!! Boundary Layer
664	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! Coefficients B for q from the Planetary
665	!! Boundary Layer
666	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
667
668
669	!! 0.2 Output variables
670
671	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: coastalflow !! Outflow on coastal points by small basins.
672	!! This is the water which flows in a disperse
673	!! way into the ocean
674	!! @tex $(kg dt_routing^{-1})$ @endtex
675	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: riverflow !! Outflow of the major rivers.
676	!! The flux will be located on the continental
677	!! grid but this should be a coastal point
678	!! @tex $(kg dt_routing^{-1})$ @endtex
679	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature
680	!! @tex $(W m^{-2})$ @endtex
681	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation
682	!! @tex $(kg m^{-2} days^{-1})$ @endtex
683
684	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf_out !! Surface specific humidity
685	!! @tex $(kg kg^{-1})$ @endtex
686	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m_out !! Surface roughness momentum (output diagnostic, m)
687	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h_out !! Surface roughness heat (output diagnostic, m)
688	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Surface albedo for visible and near-infrared
689	!! (unitless, 0-1)
690	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux
691	!! @tex $(W m^{-2})$ @endtex
692	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux
693	!! @tex $(W m^{-2})$ @endtex
694	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis_out !! Emissivity (unitless)
695	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: netco2flux !! Sum CO2 flux over PFTs
696	!! (gC/m2/dt_sechiba)
697	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_lu !! Land Cover Change CO2 flux (gC/m2/one_day)
698	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_wh !! Wood harvest CO2 flux (gC/m2/one_day)
699	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_ha !! Crop harvest CO2 flux (gC/m2/one_day)
700	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget_out !! Fraction of vegetation type (unitless, 0-1)
701	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: lai_out !! Leaf area index (m^2 m^{-2})
702	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height_out !! Vegetation Height (m)
703
704	!! 0.3 Modified
705
706	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: tq_cdrag !! Surface drag coefficient (-)
707	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! New ground temperature (K)
708
709	!! 0.4 local variables
710
711	INTEGER(i_std) :: ji, jv !! Index (unitless)
712	REAL(r_std), DIMENSION(kjpindex) :: histvar !! Computations for history files (unitless)
713	REAL(r_std), DIMENSION(kjpindex,nlut) :: histvar2 !! Computations for history files (unitless)
714	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
715	REAL(r_std), DIMENSION(kjpindex) :: sum_treefrac !! Total fraction occupied by trees (0-1, uniless)
716	REAL(r_std), DIMENSION(kjpindex) :: sum_grassfracC3 !! Total fraction occupied by C3 grasses (0-1, unitless)
717	REAL(r_std), DIMENSION(kjpindex) :: sum_grassfracC4 !! Total fraction occupied by C4 grasses (0-1, unitless)
718	REAL(r_std), DIMENSION(kjpindex) :: sum_cropfracC3 !! Total fraction occupied by C3 crops (0-1, unitess)
719	REAL(r_std), DIMENSION(kjpindex) :: sum_cropfracC4 !! Total fraction occupied by C4 crops (0-1, unitess)
720	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracNdlEvg!! Total fraction occupied by treeFracNdlEvg (0-1, unitess)
721	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracBdlEvg!! Total fraction occupied by treeFracBdlEvg (0-1, unitess)
722	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracNdlDcd!! Total fraction occupied by treeFracNdlDcd (0-1, unitess)
723	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracBdlDcd!! Total fraction occupied by treeFracBdlDcd (0-1, unitess)
724	REAL(r_std), DIMENSION(kjpindex) :: grndflux !! Net energy into soil (W/m2)
725	REAL(r_std), DIMENSION(kjpindex,nsnow) :: snowliq !! Liquid water content (m)
726	REAL(r_std), DIMENSION(kjpindex) :: snow_age_diag !! Only for diag, contains xios_default_val
727	REAL(r_std), DIMENSION(kjpindex,nnobio) :: snow_nobio_age_diag !! Only for diag, contains xios_default_val
728	REAL(r_std), DIMENSION(kjpindex) :: snowage_glob !! Snow age on total area including snow on vegetated and bare soil and nobio area @tex ($d$) @endtex
729	REAL(r_std), DIMENSION(kjpindex,nlut) :: gpplut !! GPP on landuse tile, only for diagnostics
730
731
732	!_ ================================================================================================================================
733
734	IF (printlev_loc>=3) WRITE(numout,*) 'Start sechiba_main kjpindex =',kjpindex
735
736	!! 1. Initialize variables at each time step
737	CALL sechiba_var_init (kjpindex, rau, pb, temp_air)
738
739	!! 2. Compute diffusion coefficients
740	CALL diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, &
741	& zlev, z0m, z0h, roughheight, temp_sol, temp_air, temp_growth, rau, tq_cdrag, qsurf, qair, pb , &
742	& evap_bare_lim, evap_bare_lim_ns, evapot, evapot_corr, snow, flood_frac, flood_res, &
743	& frac_nobio, snow_nobio, totfrac_nobio, &
744	& swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
745	& vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, gsmean, rveget, rstruct, cimean, gpp, &
746	& lalo, neighbours, resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
747	& hist_id, hist2_id)
748
749	!! 3. Compute energy balance
750	CALL enerbil_main (kjit, kjpindex, &
751	& index, indexveg, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, &
752	& qair, peqAcoef, peqBcoef, pb, rau, vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
753	& emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
754	& vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo, temp_sol, tsol_rad, &
755	& temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id, &
756	& precip_rain, pgflux, snowdz, temp_sol_add)
757
758
759	!! 4. Compute hydrology
760	!! 4.1 Water balance from CWRR module (11 soil layers)
761	CALL hydrol_main (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjit, kjpindex, &
762	& index, indexveg, indexsoil, indexlayer, indexnslm, &
763	& temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, &
764	& qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age, &
765	& tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, &
766	& humrel, vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, evap_bare_lim_ns, flood_frac, flood_res, &
767	& shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope_soil,&
768	& rest_id, hist_id, hist2_id,&
769	& contfrac, stempdiag, &
770	& temp_air, pb, u, v, tq_cdrag, swnet, pgflux, &
771	& snowrho,snowtemp,snowgrain,snowdz,snowheat,snowliq, &
772	& grndflux,gtemp,tot_bare_soil, &
773	& lambda_snow,cgrnd_snow,dgrnd_snow,frac_snow_veg,temp_sol_add, &
774	& mc_layh, mcl_layh, soilmoist, root_deficit)
775
776
777	!! 6. Compute surface variables (emissivity, albedo and roughness)
778	CALL condveg_main (kjit, kjpindex, index, rest_id, hist_id, hist2_id, &
779	lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
780	zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
781	drysoil_frac, height, snowdz, snowrho, tot_bare_soil, &
782	temp_air, pb, u, v, lai, &
783	emis, albedo, z0m, z0h, roughheight, &
784	frac_snow_veg, frac_snow_nobio)
785
786
787	!! 7. Compute soil thermodynamics
788	CALL thermosoil_main (kjit, kjpindex, &
789	index, indexgrnd, mcs, &
790	temp_sol_new, snow, soilcap, soilflx, &
791	shumdiag_perma, stempdiag, ftempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
792	snowdz,snowrho,snowtemp,gtemp,pb,&
793	mc_layh, mcl_layh, soilmoist, njsc,frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
794	lambda_snow, cgrnd_snow, dgrnd_snow)
795
796
797	!! 8. Compute river routing
798	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
799	!! 8.1 River routing
800	CALL routing_wrapper_main (kjit, kjpindex, index, &
801	& lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
802	& drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, stempdiag, &
803	& ftempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id, &
804	& soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw)
805	ELSE
806	!! 8.2 No routing, set variables to zero
807	riverflow(:) = zero
808	coastalflow(:) = zero
809	returnflow(:) = zero
810	reinfiltration(:) = zero
811	irrigation(:) = zero
812	flood_frac(:) = zero
813	flood_res(:) = zero
814
815	CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba)
816	CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba)
817	ENDIF
818
819	!! 9. Compute slow processes (i.e. 'daily' and annual time step)
820	CALL slowproc_main (kjit, kjpij, kjpindex, &
821	index, indexveg, lalo, neighbours, resolution, contfrac, soiltile, fraclut, nwdFraclut, &
822	temp_air, temp_sol, stempdiag, &
823	vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
824	deadleaf_cover, &
825	assim_param, &
826	lai, frac_age, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
827	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
828	co2_flux, fco2_lu, fco2_wh, fco2_ha, temp_growth, tot_bare_soil, &
829	irrigated_next, irrig_frac, reinf_slope, reinf_slope_soil)
830
831
832	!! 9.2 Compute global CO2 flux
833	netco2flux(:) = zero
834	DO jv = 2,nvm
835	netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*(1-totfrac_nobio)
836	ENDDO
837
838	!! 10. Update the temperature (temp_sol) with newly computed values
839	CALL sechiba_end (kjpindex, temp_sol_new, temp_sol)
840
841
842	!! 11. Write internal variables to output fields
843	z0m_out(:) = z0m(:)
844	z0h_out(:) = z0h(:)
845	emis_out(:) = emis(:)
846	qsurf_out(:) = qsurf(:)
847	veget_out(:,:) = veget(:,:)
848	lai_out(:,:) = lai(:,:)
849	height_out(:,:) = height(:,:)
850
851	!! 12. Write global variables to history files
852	sum_treefrac(:) = zero
853	sum_grassfracC3(:) = zero
854	sum_grassfracC4(:) = zero
855	sum_cropfracC3(:) = zero
856	sum_cropfracC4(:) = zero
857	sum_treeFracNdlEvg(:) = zero
858	sum_treeFracBdlEvg(:) = zero
859	sum_treeFracNdlDcd(:) = zero
860	sum_treeFracBdlDcd(:) = zero
861
862	DO jv = 2, nvm
863	IF (is_tree(jv) .AND. natural(jv)) THEN
864	sum_treefrac(:) = sum_treefrac(:) + veget_max(:,jv)
865	ELSE IF ((.NOT. is_tree(jv)) .AND. natural(jv)) THEN
866	! Grass
867	IF (is_c4(jv)) THEN
868	sum_grassfracC4(:) = sum_grassfracC4(:) + veget_max(:,jv)
869	ELSE
870	sum_grassfracC3(:) = sum_grassfracC3(:) + veget_max(:,jv)
871	END IF
872	ELSE
873	! Crop and trees not natural
874	IF (is_c4(jv)) THEN
875	sum_cropfracC4(:) = sum_cropfracC4(:) + veget_max(:,jv)
876	ELSE
877	sum_cropfracC3(:) = sum_cropfracC3(:) + veget_max(:,jv)
878	END IF
879	ENDIF
880
881	IF (is_tree(jv)) THEN
882	IF (is_evergreen(jv)) THEN
883	IF (is_needleleaf(jv)) THEN
884	! Fraction for needleleaf evergreen trees (treeFracNdlEvg)
885	sum_treeFracNdlEvg(:) = sum_treeFracNdlEvg(:) + veget_max(:,jv)
886	ELSE
887	! Fraction for broadleaf evergreen trees (treeFracBdlEvg)
888	sum_treeFracBdlEvg(:) = sum_treeFracBdlEvg(:) + veget_max(:,jv)
889	END IF
890	ELSE IF (is_deciduous(jv)) THEN
891	IF (is_needleleaf(jv)) THEN
892	! Fraction for needleleaf deciduous trees (treeFracNdlDcd)
893	sum_treeFracNdlDcd(:) = sum_treeFracNdlDcd(:) + veget_max(:,jv)
894	ELSE
895	! Fraction for broadleafs deciduous trees (treeFracBdlDcd)
896	sum_treeFracBdlDcd(:) = sum_treeFracBdlDcd(:) + veget_max(:,jv)
897	END IF
898	END IF
899	END IF
900	ENDDO
901
902	histvar(:)=zero
903	DO jv = 2, nvm
904	IF (is_deciduous(jv)) THEN
905	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
906	ENDIF
907	ENDDO
908	CALL xios_orchidee_send_field("treeFracPrimDec",histvar)
909
910	histvar(:)=zero
911	DO jv = 2, nvm
912	IF (is_evergreen(jv)) THEN
913	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
914	ENDIF
915	ENDDO
916	CALL xios_orchidee_send_field("treeFracPrimEver",histvar)
917
918	histvar(:)=zero
919	DO jv = 2, nvm
920	IF ( .NOT.(is_c4(jv)) ) THEN
921	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
922	ENDIF
923	ENDDO
924	CALL xios_orchidee_send_field("c3PftFrac",histvar)
925
926	histvar(:)=zero
927	DO jv = 2, nvm
928	IF ( is_c4(jv) ) THEN
929	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
930	ENDIF
931	ENDDO
932	CALL xios_orchidee_send_field("c4PftFrac",histvar)
933
934	CALL xios_orchidee_send_field("temp_sol_new",temp_sol_new)
935	CALL xios_orchidee_send_field("fluxsens",fluxsens)
936	CALL xios_orchidee_send_field("fluxlat",fluxlat)
937
938
939	! Add XIOS default value where no snow
940	DO ji=1,kjpindex
941	IF (snow(ji) .GT. zero) THEN
942	snow_age_diag(ji) = snow_age(ji)
943	snow_nobio_age_diag(ji,:) = snow_nobio_age(ji,:)
944
945	snowage_glob(ji) = snow_age(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji)) + &
946	SUM(snow_nobio_age(ji,:)frac_snow_nobio(ji,:)frac_nobio(ji,:))
947	IF (snowage_glob(ji) .NE. 0) snowage_glob(ji) = snowage_glob(ji) / &
948	(frac_snow_veg(ji)(1-totfrac_nobio(ji)) + SUM(frac_snow_nobio(ji,:)frac_nobio(ji,:)))
949	ELSE
950	snow_age_diag(ji) = xios_default_val
951	snow_nobio_age_diag(ji,:) = xios_default_val
952	snowage_glob(ji) = xios_default_val
953	END IF
954	END DO
955
956	CALL xios_orchidee_send_field("snow",snow)
957	CALL xios_orchidee_send_field("snowage",snow_age_diag)
958	CALL xios_orchidee_send_field("snownobio",snow_nobio)
959	CALL xios_orchidee_send_field("snownobioage",snow_nobio_age_diag)
960	CALL xios_orchidee_send_field("snowage_glob",snowage_glob)
961
962	CALL xios_orchidee_send_field("frac_snow", SUM(frac_snow_nobio,2)totfrac_nobio+frac_snow_veg(1-totfrac_nobio))
963	CALL xios_orchidee_send_field("frac_snow_veg", frac_snow_veg)
964	CALL xios_orchidee_send_field("frac_snow_nobio", frac_snow_nobio)
965	CALL xios_orchidee_send_field("reinf_slope",reinf_slope)
966	CALL xios_orchidee_send_field("njsc",REAL(njsc, r_std))
967	CALL xios_orchidee_send_field("vegetfrac",veget)
968	CALL xios_orchidee_send_field("maxvegetfrac",veget_max)
969	CALL xios_orchidee_send_field("irrigmap_dyn",irrigated_next)
970	CALL xios_orchidee_send_field("aei_sw",fraction_aeirrig_sw)
971	CALL xios_orchidee_send_field("nobiofrac",frac_nobio)
972	CALL xios_orchidee_send_field("soiltile",soiltile)
973	CALL xios_orchidee_send_field("rstruct",rstruct)
974	CALL xios_orchidee_send_field("gpp",gpp/dt_sechiba)
975	CALL xios_orchidee_send_field("gpp_ipcc2",SUM(gpp,dim=2)/dt_sechiba)
976
977	histvar(:)=zero
978	DO jv = 2, nvm
979	IF ( .NOT. is_tree(jv) .AND. natural(jv) ) THEN
980	histvar(:) = histvar(:) + gpp(:,jv)
981	ENDIF
982	ENDDO
983	CALL xios_orchidee_send_field("gppgrass",histvar/dt_sechiba)
984
985	histvar(:)=zero
986	DO jv = 2, nvm
987	IF ( (.NOT. is_tree(jv)) .AND. (.NOT. natural(jv)) ) THEN
988	histvar(:) = histvar(:) + gpp(:,jv)
989	ENDIF
990	ENDDO
991	CALL xios_orchidee_send_field("gppcrop",histvar/dt_sechiba)
992
993	histvar(:)=zero
994	DO jv = 2, nvm
995	IF ( is_tree(jv) ) THEN
996	histvar(:) = histvar(:) + gpp(:,jv)
997	ENDIF
998	ENDDO
999	CALL xios_orchidee_send_field("gpptree",histvar/dt_sechiba)
1000	CALL xios_orchidee_send_field("nee",co2_flux/1.e3/one_day)
1001	CALL xios_orchidee_send_field("drysoil_frac",drysoil_frac)
1002	CALL xios_orchidee_send_field("vevapflo",vevapflo/dt_sechiba)
1003	CALL xios_orchidee_send_field("k_litt",k_litt)
1004	CALL xios_orchidee_send_field("beta",vbeta)
1005	CALL xios_orchidee_send_field("vbeta1",vbeta1)
1006	CALL xios_orchidee_send_field("vbeta2",vbeta2)
1007	CALL xios_orchidee_send_field("vbeta3",vbeta3)
1008	CALL xios_orchidee_send_field("vbeta4",vbeta4)
1009	CALL xios_orchidee_send_field("vbeta5",vbeta5)
1010	CALL xios_orchidee_send_field("gsmean",gsmean)
1011	CALL xios_orchidee_send_field("cimean",cimean)
1012	CALL xios_orchidee_send_field("rveget",rveget)
1013
1014	histvar(:)=SUM(vevapwet(:,:),dim=2)
1015	CALL xios_orchidee_send_field("evspsblveg",histvar/dt_sechiba)
1016	histvar(:)= vevapnu(:)+vevapsno(:)
1017	CALL xios_orchidee_send_field("evspsblsoi",histvar/dt_sechiba)
1018	histvar(:)=SUM(transpir(:,:),dim=2)
1019	CALL xios_orchidee_send_field("tran",histvar/dt_sechiba)
1020
1021	! For CMIP6 data request: the following fractions are fractions of the total grid-cell,
1022	! which explains the multiplication by contfrac
1023	CALL xios_orchidee_send_field("treeFrac",sum_treefrac(:)100contfrac(:))
1024	CALL xios_orchidee_send_field("grassFracC3",sum_grassfracC3(:)100contfrac(:))
1025	CALL xios_orchidee_send_field("grassFracC4",sum_grassfracC4(:)100contfrac(:))
1026	CALL xios_orchidee_send_field("cropFracC3",sum_cropfracC3(:)100contfrac(:))
1027	CALL xios_orchidee_send_field("cropFracC4",sum_cropfracC4(:)100contfrac(:))
1028	CALL xios_orchidee_send_field("treeFracNdlEvg",sum_treeFracNdlEvg(:)100contfrac(:))
1029	CALL xios_orchidee_send_field("treeFracBdlEvg",sum_treeFracBdlEvg(:)100contfrac(:))
1030	CALL xios_orchidee_send_field("treeFracNdlDcd",sum_treeFracNdlDcd(:)100contfrac(:))
1031	CALL xios_orchidee_send_field("treeFracBdlDcd",sum_treeFracBdlDcd(:)100contfrac(:))
1032
1033	histvar(:)=veget_max(:,1)100contfrac(:)
1034	CALL xios_orchidee_send_field("baresoilFrac",histvar)
1035	histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)100contfrac(:)
1036	CALL xios_orchidee_send_field("residualFrac",histvar)
1037
1038	! For CMIP6 data request: cnc = canopy cover fraction over land area
1039	histvar(:)=zero
1040	DO jv=2,nvm
1041	histvar(:) = histvar(:) + veget_max(:,jv)*100
1042	END DO
1043	CALL xios_orchidee_send_field("cnc",histvar)
1044
1045	CALL xios_orchidee_send_field("tsol_rad",tsol_rad-273.15)
1046	CALL xios_orchidee_send_field("qsurf",qsurf)
1047	CALL xios_orchidee_send_field("emis",emis)
1048	CALL xios_orchidee_send_field("z0m",z0m)
1049	CALL xios_orchidee_send_field("z0h",z0h)
1050	CALL xios_orchidee_send_field("roughheight",roughheight)
1051	CALL xios_orchidee_send_field("lai",lai)
1052	histvar(:)=zero
1053	DO ji = 1, kjpindex
1054	IF (SUM(veget_max(ji,:)) > zero) THEN
1055	DO jv=2,nvm
1056	histvar(ji) = histvar(ji) + veget_max(ji,jv)*lai(ji,jv)/SUM(veget_max(ji,:))
1057	END DO
1058	END IF
1059	END DO
1060	CALL xios_orchidee_send_field("LAImean",histvar)
1061	CALL xios_orchidee_send_field("vevapsno",vevapsno/dt_sechiba)
1062	CALL xios_orchidee_send_field("vevapp",vevapp/dt_sechiba)
1063	CALL xios_orchidee_send_field("vevapnu",vevapnu/dt_sechiba)
1064	CALL xios_orchidee_send_field("transpir",transpir*one_day/dt_sechiba)
1065	CALL xios_orchidee_send_field("transpot",transpot*one_day/dt_sechiba)
1066	CALL xios_orchidee_send_field("inter",vevapwet*one_day/dt_sechiba)
1067	histvar(:)=zero
1068	DO jv=1,nvm
1069	histvar(:) = histvar(:) + vevapwet(:,jv)
1070	ENDDO
1071	CALL xios_orchidee_send_field("ECanop",histvar/dt_sechiba)
1072	histvar(:)=zero
1073	DO jv=1,nvm
1074	histvar(:) = histvar(:) + transpir(:,jv)
1075	ENDDO
1076	CALL xios_orchidee_send_field("TVeg",histvar/dt_sechiba)
1077
1078
1079	!! Calculate diagnostic variables on Landuse tiles for LUMIP/CMIP6
1080
1081	! Calculate fraction of landuse tiles related to the whole grid cell
1082	DO jv=1,nlut
1083	histvar2(:,jv) = fraclut(:,jv) * contfrac(:)
1084	END DO
1085	CALL xios_orchidee_send_field("fraclut",histvar2)
1086
1087	CALL xios_orchidee_send_field("nwdFraclut",nwdFraclut(:,:))
1088
1089	! Calculate GPP on landuse tiles
1090	! val_exp is used as missing value where the values are not known i.e. where the tile is not represented
1091	! or for pasture (id_pst) or urban land (id_urb).
1092	gpplut(:,:)=0
1093	DO jv=1,nvm
1094	IF (natural(jv)) THEN
1095	gpplut(:,id_psl) = gpplut(:,id_psl) + gpp(:,jv)
1096	ELSE
1097	gpplut(:,id_crp) = gpplut(:,id_crp) + gpp(:,jv)
1098	ENDIF
1099	END DO
1100
1101	! Transform from gC/m2/s into kgC/m2/s
1102	WHERE (fraclut(:,id_psl)>min_sechiba)
1103	gpplut(:,id_psl) = gpplut(:,id_psl)/fraclut(:,id_psl)/1000
1104	ELSEWHERE
1105	gpplut(:,id_psl) = xios_default_val
1106	END WHERE
1107	WHERE (fraclut(:,id_crp)>min_sechiba)
1108	gpplut(:,id_crp) = gpplut(:,id_crp)/fraclut(:,id_crp)/1000
1109	ELSEWHERE
1110	gpplut(:,id_crp) = xios_default_val
1111	END WHERE
1112	gpplut(:,id_pst) = xios_default_val
1113	gpplut(:,id_urb) = xios_default_val
1114
1115	CALL xios_orchidee_send_field("gpplut",gpplut)
1116
1117
1118	IF ( .NOT. almaoutput ) THEN
1119	! Write history file in IPSL-format
1120	CALL histwrite_p(hist_id, 'beta', kjit, vbeta, kjpindex, index)
1121	CALL histwrite_p(hist_id, 'z0m', kjit, z0m, kjpindex, index)
1122	CALL histwrite_p(hist_id, 'z0h', kjit, z0h, kjpindex, index)
1123	CALL histwrite_p(hist_id, 'roughheight', kjit, roughheight, kjpindex, index)
1124	CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1125	CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1126	CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1127	CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1128	CALL histwrite_p(hist_id, 'subli', kjit, vevapsno, kjpindex, index)
1129	CALL histwrite_p(hist_id, 'evapnu', kjit, vevapnu, kjpindex, index)
1130	CALL histwrite_p(hist_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
1131	CALL histwrite_p(hist_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
1132	CALL histwrite_p(hist_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
1133	CALL histwrite_p(hist_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
1134	CALL histwrite_p(hist_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
1135	CALL histwrite_p(hist_id, 'vbeta4', kjit, vbeta4, kjpindex, index)
1136	CALL histwrite_p(hist_id, 'vbeta5', kjit, vbeta5, kjpindex, index)
1137	CALL histwrite_p(hist_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
1138	CALL histwrite_p(hist_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
1139	CALL histwrite_p(hist_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
1140	CALL histwrite_p(hist_id, 'snow', kjit, snow, kjpindex, index)
1141	CALL histwrite_p(hist_id, 'snowage', kjit, snow_age, kjpindex, index)
1142	CALL histwrite_p(hist_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
1143	CALL histwrite_p(hist_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
1144
1145	CALL histwrite_p(hist_id, 'grndflux', kjit, grndflux, kjpindex,index)
1146	CALL histwrite_p(hist_id, 'snowtemp',kjit,snowtemp,kjpindex*nsnow,indexsnow)
1147	CALL histwrite_p(hist_id, 'snowliq', kjit,snowliq,kjpindex*nsnow,indexsnow)
1148	CALL histwrite_p(hist_id, 'snowdz', kjit,snowdz,kjpindex*nsnow,indexsnow)
1149	CALL histwrite_p(hist_id, 'snowrho', kjit,snowrho,kjpindex*nsnow,indexsnow)
1150	CALL histwrite_p(hist_id, 'snowgrain',kjit,snowgrain,kjpindex*nsnow,indexsnow)
1151	CALL histwrite_p(hist_id, 'snowheat',kjit,snowheat,kjpindex*nsnow,indexsnow)
1152
1153	CALL histwrite_p(hist_id,'frac_snow_veg',kjit,frac_snow_veg,kjpindex,index)
1154	CALL histwrite_p(hist_id, 'frac_snow_nobio', kjit,frac_snow_nobio,kjpindex*nnobio, indexnobio)
1155	CALL histwrite_p(hist_id, 'pgflux',kjit,pgflux,kjpindex,index)
1156	CALL histwrite_p(hist_id, 'soiltile', kjit, soiltile, kjpindex*nstm, indexsoil)
1157
1158	CALL histwrite_p(hist_id, 'soilindex', kjit, REAL(njsc, r_std), kjpindex, index)
1159	CALL histwrite_p(hist_id, 'reinf_slope', kjit, reinf_slope, kjpindex, index)
1160	CALL histwrite_p(hist_id, 'k_litt', kjit, k_litt, kjpindex, index)
1161
1162	IF ( do_floodplains ) THEN
1163	CALL histwrite_p(hist_id, 'evapflo', kjit, vevapflo, kjpindex, index)
1164	CALL histwrite_p(hist_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
1165	ENDIF
1166
1167	CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)
1168	CALL histwrite_p(hist_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
1169	CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
1170
1171	IF ( ok_stomate ) THEN
1172	CALL histwrite_p(hist_id, 'nee', kjit, co2_flux/1.e3/one_day, kjpindex*nvm, indexveg)
1173	ENDIF
1174
1175	histvar(:)=SUM(vevapwet(:,:),dim=2)
1176	CALL histwrite_p(hist_id, 'evspsblveg', kjit, histvar, kjpindex, index)
1177
1178	histvar(:)= vevapnu(:)+vevapsno(:)
1179	CALL histwrite_p(hist_id, 'evspsblsoi', kjit, histvar, kjpindex, index)
1180
1181	histvar(:)=SUM(transpir(:,:),dim=2)
1182	CALL histwrite_p(hist_id, 'tran', kjit, histvar, kjpindex, index)
1183
1184	histvar(:)= sum_treefrac(:)100contfrac(:)
1185	CALL histwrite_p(hist_id, 'treeFrac', kjit, histvar, kjpindex, index)
1186
1187	histvar(:)= (sum_grassfracC3(:)+sum_grassfracC4(:))100contfrac(:)
1188	CALL histwrite_p(hist_id, 'grassFrac', kjit, histvar, kjpindex, index)
1189
1190	histvar(:)= (sum_cropfracC3(:)+sum_cropfracC4(:))100contfrac(:)
1191	CALL histwrite_p(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)
1192
1193	histvar(:)=veget_max(:,1)100contfrac(:)
1194	CALL histwrite_p(hist_id, 'baresoilFrac', kjit, histvar, kjpindex, index)
1195
1196	histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)100contfrac(:)
1197	CALL histwrite_p(hist_id, 'residualFrac', kjit, histvar, kjpindex, index)
1198	ELSE
1199	! Write history file in ALMA format
1200	CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1201	CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1202	CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1203	CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1204	CALL histwrite_p(hist_id, 'ESoil', kjit, vevapnu, kjpindex, index)
1205	CALL histwrite_p(hist_id, 'EWater', kjit, vevapflo, kjpindex, index)
1206	CALL histwrite_p(hist_id, 'SWE', kjit, snow, kjpindex, index)
1207	histvar(:)=zero
1208	DO jv=1,nvm
1209	histvar(:) = histvar(:) + transpir(:,jv)
1210	ENDDO
1211	CALL histwrite_p(hist_id, 'TVeg', kjit, histvar, kjpindex, index)
1212	histvar(:)=zero
1213	DO jv=1,nvm
1214	histvar(:) = histvar(:) + vevapwet(:,jv)
1215	ENDDO
1216	CALL histwrite_p(hist_id, 'ECanop', kjit, histvar, kjpindex, index)
1217	CALL histwrite_p(hist_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
1218	!
1219	CALL histwrite_p(hist_id, 'Z0m', kjit, z0m, kjpindex, index)
1220	CALL histwrite_p(hist_id, 'Z0h', kjit, z0h, kjpindex, index)
1221	CALL histwrite_p(hist_id, 'EffectHeight', kjit, roughheight, kjpindex, index)
1222	!
1223	IF ( do_floodplains ) THEN
1224	CALL histwrite_p(hist_id, 'Qflood', kjit, vevapflo, kjpindex, index)
1225	CALL histwrite_p(hist_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
1226	ENDIF
1227
1228	CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)
1229	CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
1230	CALL histwrite_p(hist_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
1231
1232	IF ( ok_stomate ) THEN
1233	CALL histwrite_p(hist_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)
1234	ENDIF
1235	ENDIF ! almaoutput
1236
1237	!! 13. Write additional output file with higher frequency
1238	IF ( hist2_id > 0 ) THEN
1239	IF ( .NOT. almaoutput ) THEN
1240	! Write history file in IPSL-format
1241	CALL histwrite_p(hist2_id, 'tsol_rad', kjit, tsol_rad, kjpindex, index)
1242	CALL histwrite_p(hist2_id, 'qsurf', kjit, qsurf, kjpindex, index)
1243	CALL histwrite_p(hist2_id, 'albedo', kjit, albedo, kjpindex*2, indexalb)
1244	CALL histwrite_p(hist2_id, 'emis', kjit, emis, kjpindex, index)
1245	CALL histwrite_p(hist2_id, 'beta', kjit, vbeta, kjpindex, index)
1246	CALL histwrite_p(hist2_id, 'z0m', kjit, z0m, kjpindex, index)
1247	CALL histwrite_p(hist2_id, 'z0h', kjit, z0h, kjpindex, index)
1248	CALL histwrite_p(hist2_id, 'roughheight', kjit, roughheight, kjpindex, index)
1249	CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1250	CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1251	CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1252	CALL histwrite_p(hist2_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1253	CALL histwrite_p(hist2_id, 'subli', kjit, vevapsno, kjpindex, index)
1254	IF ( do_floodplains ) THEN
1255	CALL histwrite_p(hist2_id, 'vevapflo', kjit, vevapflo, kjpindex, index)
1256	CALL histwrite_p(hist2_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
1257	ENDIF
1258	CALL histwrite_p(hist2_id, 'vevapnu', kjit, vevapnu, kjpindex, index)
1259	CALL histwrite_p(hist2_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
1260	CALL histwrite_p(hist2_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
1261	CALL histwrite_p(hist2_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
1262	CALL histwrite_p(hist2_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
1263	CALL histwrite_p(hist2_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
1264	CALL histwrite_p(hist2_id, 'vbeta4', kjit, vbeta4, kjpindex, index)
1265	CALL histwrite_p(hist2_id, 'vbeta5', kjit, vbeta5, kjpindex, index)
1266	CALL histwrite_p(hist2_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
1267	CALL histwrite_p(hist2_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
1268	CALL histwrite_p(hist2_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
1269	CALL histwrite_p(hist2_id, 'snow', kjit, snow, kjpindex, index)
1270	CALL histwrite_p(hist2_id, 'snowage', kjit, snow_age, kjpindex, index)
1271	CALL histwrite_p(hist2_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
1272	CALL histwrite_p(hist2_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
1273	CALL histwrite_p(hist2_id, 'soilindex', kjit, REAL(njsc, r_std), kjpindex, index)
1274	CALL histwrite_p(hist2_id, 'reinf_slope', kjit, reinf_slope, kjpindex, index)
1275
1276	CALL histwrite_p(hist2_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)
1277	CALL histwrite_p(hist2_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
1278	CALL histwrite_p(hist2_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
1279
1280	IF ( ok_stomate ) THEN
1281	CALL histwrite_p(hist2_id, 'nee', kjit, co2_flux/1.e3/one_day, kjpindex*nvm, indexveg)
1282	ENDIF
1283	ELSE
1284	! Write history file in ALMA format
1285	CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1286	CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1287	CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1288	CALL histwrite_p(hist2_id, 'ESoil', kjit, vevapnu, kjpindex, index)
1289	IF ( do_floodplains ) THEN
1290	CALL histwrite_p(hist2_id, 'EWater', kjit, vevapflo, kjpindex, index)
1291	CALL histwrite_p(hist2_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
1292	ENDIF
1293	CALL histwrite_p(hist2_id, 'SWE', kjit, snow, kjpindex, index)
1294	histvar(:)=zero
1295	DO jv=1,nvm
1296	histvar(:) = histvar(:) + transpir(:,jv)
1297	ENDDO
1298	CALL histwrite_p(hist2_id, 'TVeg', kjit, histvar, kjpindex, index)
1299	histvar(:)=zero
1300	DO jv=1,nvm
1301	histvar(:) = histvar(:) + vevapwet(:,jv)
1302	ENDDO
1303	CALL histwrite_p(hist2_id, 'ECanop', kjit, histvar, kjpindex, index)
1304	CALL histwrite_p(hist2_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
1305	CALL histwrite_p(hist2_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
1306
1307	IF ( ok_stomate ) THEN
1308	CALL histwrite_p(hist2_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)
1309	ENDIF
1310	ENDIF ! almaoutput
1311	ENDIF ! hist2_id
1312
1313
1314	!! Change the vegetation fractions if a new map was read in slowproc. This is done
1315	!! after lcchange has been done in stomatelpj
1316	IF (done_stomate_lcchange) THEN
1317	CALL slowproc_change_frac(kjpindex, lai, &
1318	veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile, fraclut, nwdFraclut)
1319	done_stomate_lcchange=.FALSE.
1320	END IF
1321
1322	!! 14. If it is the last time step, write restart files
1323	IF (ldrestart_write) THEN
1324	CALL sechiba_finalize( &
1325	kjit, kjpij, kjpindex, index, rest_id, &
1326	tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad)
1327	END IF
1328
1329	END SUBROUTINE sechiba_main
1330
1331
1332	!! =============================================================================================================================
1333	!! SUBROUTINE: sechiba_finalize
1334	!!
1335	!>\BRIEF Finalize all modules by calling their "_finalize" subroutines.
1336	!!
1337	!! DESCRIPTION: Finalize all modules by calling their "_finalize" subroutines. These subroutines will write variables to
1338	!! restart file.
1339	!!
1340	!! \n
1341	!_ ==============================================================================================================================
1342
1343	SUBROUTINE sechiba_finalize( &
1344	kjit, kjpij, kjpindex, index, rest_id, &
1345	tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad)
1346
1347	!! 0.1 Input variables
1348	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
1349	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
1350	!! (unitless)
1351	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
1352	!! (unitless)
1353	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
1354	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the pixels on the map.
1355	!! Sechiba uses a reduced grid excluding oceans
1356	!! ::index contains the indices of the
1357	!! terrestrial pixels only! (unitless)
1358	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tsol_rad !! Radiative surface temperature
1359	!! @tex $(W m^{-2})$ @endtex
1360	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapp !! Total of evaporation
1361	!! @tex $(kg m^{-2} days^{-1})$ @endtex
1362	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxsens !! Sensible heat flux
1363	!! @tex $(W m^{-2})$ @endtex
1364	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxlat !! Latent heat flux
1365	!! @tex $(W m^{-2})$ @endtex
1366	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tq_cdrag !! Surface drag coefficient (-)
1367
1368	!! 0.2 Local variables
1369	INTEGER(i_std) :: ji, jv !! Index (unitless)
1370	REAL(r_std), DIMENSION(kjpindex) :: histvar !! Computations for history files (unitless)
1371	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
1372
1373
1374	!! Write restart file for the next simulation from SECHIBA and other modules
1375
1376	IF (printlev_loc>=3) WRITE (numout,*) 'Start sechiba_finalize for writing restart files'
1377
1378	!! 1. Call diffuco_finalize to write restart files
1379	CALL diffuco_finalize (kjit, kjpindex, rest_id, rstruct )
1380
1381	!! 2. Call energy budget to write restart files
1382	CALL enerbil_finalize (kjit, kjpindex, rest_id, &
1383	evapot, evapot_corr, temp_sol, tsol_rad, &
1384	qsurf, fluxsens, fluxlat, vevapp )
1385
1386	!! 3. Call hydrology to write restart files
1387	CALL hydrol_finalize( kjit, kjpindex, rest_id, vegstress, &
1388	qsintveg, humrel, snow, snow_age, snow_nobio, &
1389	snow_nobio_age, snowrho, snowtemp, snowdz, &
1390	snowheat, snowgrain, &
1391	drysoil_frac, evap_bare_lim, evap_bare_lim_ns)
1392
1393	!! 4. Call condveg to write surface variables to restart files
1394	CALL condveg_finalize (kjit, kjpindex, rest_id, z0m, z0h, roughheight)
1395
1396	!! 5. Call soil thermodynamic to write restart files
1397	CALL thermosoil_finalize (kjit, kjpindex, rest_id, gtemp, &
1398	soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
1399
1400
1401	!! 6. Add river routing to restart files
1402	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
1403	!! 6.1 Call river routing to write restart files
1404	CALL routing_wrapper_finalize( kjit, kjpindex, rest_id, flood_frac, flood_res )
1405	ELSE
1406	!! 6.2 No routing, set variables to zero
1407	reinfiltration(:) = zero
1408	returnflow(:) = zero
1409	irrigation(:) = zero
1410	flood_frac(:) = zero
1411	flood_res(:) = zero
1412	ENDIF
1413
1414	!! 7. Call slowproc_main to add 'daily' and annual variables to restart file
1415	CALL slowproc_finalize (kjit, kjpindex, rest_id, index, &
1416	njsc, lai, height, veget, &
1417	frac_nobio, veget_max, reinf_slope, &
1418	ks, nvan, avan, mcr, &
1419	mcs, mcfc, mcw, &
1420	assim_param, frac_age, fraction_aeirrig_sw)
1421
1422	IF (printlev_loc>=3) WRITE (numout,*) 'sechiba_finalize done'
1423
1424	END SUBROUTINE sechiba_finalize
1425
1426
1427	!! ==============================================================================================================================\n
1428	!! SUBROUTINE : sechiba_init
1429	!!
1430	!>\BRIEF Dynamic allocation of the variables, the dimensions of the
1431	!! variables are determined by user-specified settings.
1432	!!
1433	!! DESCRIPTION : The domain size (:: kjpindex) is used to allocate the correct
1434	!! dimensions to all variables in sechiba. Depending on the variable, its
1435	!! dimensions are also determined by the number of PFT's (::nvm), number of
1436	!! soil types (::nstm), number of non-vegetative surface types (::nnobio),
1437	!! number of soil levels (::ngrnd), number of soil layers in the hydrological
1438	!! model (i.e. cwrr) (::nslm). Values for these variables are set in
1439	!! constantes_soil.f90 and constantes_veg.f90.\n
1440	!!
1441	!! Memory is allocated for all Sechiba variables and new indexing tables
1442	!! are build making use of both (::kjpij) and (::kjpindex). New indexing tables
1443	!! are needed because a single pixel can contain several PFTs, soil types, etc.
1444	!! The new indexing tables have separate indices for the different
1445	!! PFTs, soil types, etc.\n
1446	!!
1447	!! RECENT CHANGE(S): None
1448	!!
1449	!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
1450	!! variables. However, the routine allocates memory and builds new indexing
1451	!! variables for later use.
1452	!!
1453	!! REFERENCE(S) : None
1454	!!
1455	!! FLOWCHART : None
1456	!! \n
1457	!_ ================================================================================================================================
1458
1459	SUBROUTINE sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)
1460
1461	!! 0.1 Input variables
1462
1463	INTEGER(i_std), INTENT (in) :: kjit !! Time step number (unitless)
1464	INTEGER(i_std), INTENT (in) :: kjpij !! Total size of the un-compressed grid (unitless)
1465	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
1466	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
1467	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless)
1468	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates (latitude,longitude)
1469	!! for pixels (degrees)
1470	!! 0.2 Output variables
1471
1472	!! 0.3 Modified variables
1473
1474	!! 0.4 Local variables
1475
1476	INTEGER(i_std) :: ier !! Check errors in memory allocation (unitless)
1477	INTEGER(i_std) :: ji, jv !! Indeces (unitless)
1478	!_ ==============================================================================================================================
1479
1480	!! 1. Initialize variables
1481
1482	! Dynamic allocation with user-specified dimensions on first call
1483	IF (l_first_sechiba) THEN
1484	l_first_sechiba=.FALSE.
1485	ELSE
1486	CALL ipslerr_p(3,'sechiba_init',' l_first_sechiba false . we stop ','','')
1487	ENDIF
1488
1489	!! Initialize local printlev
1490	printlev_loc=get_printlev('sechiba')
1491
1492
1493	!! 1.1 Initialize 3D vegetation indexation table
1494	ALLOCATE (indexveg(kjpindex*nvm),stat=ier)
1495	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexveg','','')
1496
1497	ALLOCATE (indexlai(kjpindex*(nlai+1)),stat=ier)
1498	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlai','','')
1499
1500	ALLOCATE (indexsoil(kjpindex*nstm),stat=ier)
1501	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsoil','','')
1502
1503	ALLOCATE (indexnobio(kjpindex*nnobio),stat=ier)
1504	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnobio','','')
1505
1506	ALLOCATE (indexgrnd(kjpindex*ngrnd),stat=ier)
1507	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexgrnd','','')
1508
1509	ALLOCATE (indexsnow(kjpindex*nsnow),stat=ier)
1510	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsnow','','')
1511
1512	ALLOCATE (indexlayer(kjpindex*nslm),stat=ier)
1513	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlayer','','')
1514
1515	ALLOCATE (indexnslm(kjpindex*nslm),stat=ier)
1516	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnslm','','')
1517
1518	ALLOCATE (indexalb(kjpindex*2),stat=ier)
1519	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexalb','','')
1520
1521	!! 1.2 Initialize 1D array allocation with restartable value
1522	ALLOCATE (flood_res(kjpindex),stat=ier)
1523	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for flood_res','','')
1524	flood_res(:) = undef_sechiba
1525
1526	ALLOCATE (flood_frac(kjpindex),stat=ier)
1527	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for kjpindex','','')
1528	flood_frac(:) = undef_sechiba
1529
1530	ALLOCATE (snow(kjpindex),stat=ier)
1531	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow','','')
1532	snow(:) = undef_sechiba
1533
1534	ALLOCATE (snow_age(kjpindex),stat=ier)
1535	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_age','','')
1536	snow_age(:) = undef_sechiba
1537
1538	ALLOCATE (drysoil_frac(kjpindex),stat=ier)
1539	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drysoil_frac','','')
1540
1541	ALLOCATE (evap_bare_lim(kjpindex),stat=ier)
1542	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim','','')
1543
1544	ALLOCATE (evap_bare_lim_ns(kjpindex,nstm),stat=ier)
1545	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim_ns','','')
1546
1547	ALLOCATE (evapot(kjpindex),stat=ier)
1548	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot','','')
1549	evapot(:) = undef_sechiba
1550
1551	ALLOCATE (evapot_corr(kjpindex),stat=ier)
1552	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot_corr','','')
1553
1554	ALLOCATE (humrel(kjpindex,nvm),stat=ier)
1555	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for humrel','','')
1556	humrel(:,:) = undef_sechiba
1557
1558	ALLOCATE (vegstress(kjpindex,nvm),stat=ier)
1559	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vegstress','','')
1560	vegstress(:,:) = undef_sechiba
1561
1562	ALLOCATE (njsc(kjpindex),stat=ier)
1563	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for njsc','','')
1564	njsc(:)= undef_int
1565
1566	ALLOCATE (soiltile(kjpindex,nstm),stat=ier)
1567	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soiltile','','')
1568
1569	ALLOCATE (fraclut(kjpindex,nlut),stat=ier)
1570	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fraclut','','')
1571
1572	ALLOCATE (nwdFraclut(kjpindex,nlut),stat=ier)
1573	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for nwdFraclut','','')
1574
1575	ALLOCATE (reinf_slope(kjpindex),stat=ier)
1576	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope','','')
1577
1578	ALLOCATE (reinf_slope_soil(kjpindex, nstm),stat=ier)
1579	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope_soil','','') !
1580
1581	!salma: adding soil params
1582	ALLOCATE (ks(kjpindex),stat=ier)
1583	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ks','','')
1584
1585	ALLOCATE (nvan(kjpindex),stat=ier)
1586	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for nvan ','','')
1587
1588	ALLOCATE (avan(kjpindex),stat=ier)
1589	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for avan','','')
1590
1591	ALLOCATE (mcr(kjpindex),stat=ier)
1592	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcr','','')
1593
1594	ALLOCATE (mcs(kjpindex),stat=ier)
1595	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcs','','')
1596
1597	ALLOCATE (mcfc(kjpindex),stat=ier)
1598	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcfc','','')
1599
1600	ALLOCATE (mcw(kjpindex),stat=ier)
1601	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcw','','')
1602	!end salma: adding soil params
1603
1604
1605
1606
1607	ALLOCATE (vbeta1(kjpindex),stat=ier)
1608	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta1','','')
1609
1610	ALLOCATE (vbeta4(kjpindex),stat=ier)
1611	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta4','','')
1612
1613	ALLOCATE (vbeta5(kjpindex),stat=ier)
1614	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta5','','')
1615
1616	ALLOCATE (soilcap(kjpindex),stat=ier)
1617	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilcap','','')
1618
1619	ALLOCATE (soilflx(kjpindex),stat=ier)
1620	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilflx','','')
1621
1622	ALLOCATE (temp_sol(kjpindex),stat=ier)
1623	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol','','')
1624	temp_sol(:) = undef_sechiba
1625
1626	ALLOCATE (qsurf(kjpindex),stat=ier)
1627	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsurf','','')
1628	qsurf(:) = undef_sechiba
1629
1630	!! 1.3 Initialize 2D array allocation with restartable value
1631	ALLOCATE (qsintveg(kjpindex,nvm),stat=ier)
1632	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintveg','','')
1633	qsintveg(:,:) = undef_sechiba
1634
1635	ALLOCATE (vbeta2(kjpindex,nvm),stat=ier)
1636	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta2','','')
1637
1638	ALLOCATE (vbeta3(kjpindex,nvm),stat=ier)
1639	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3','','')
1640
1641	ALLOCATE (vbeta3pot(kjpindex,nvm),stat=ier)
1642	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3pot','','')
1643
1644	ALLOCATE (gsmean(kjpindex,nvm),stat=ier)
1645	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gsmean','','')
1646
1647	ALLOCATE (cimean(kjpindex,nvm),stat=ier)
1648	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cimean','','')
1649
1650	ALLOCATE (gpp(kjpindex,nvm),stat=ier)
1651	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gpp','','')
1652	gpp(:,:) = undef_sechiba
1653
1654
1655	ALLOCATE (temp_growth(kjpindex),stat=ier)
1656	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_growth','','')
1657	temp_growth(:) = undef_sechiba
1658
1659	ALLOCATE (veget(kjpindex,nvm),stat=ier)
1660	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget','','')
1661	veget(:,:)=undef_sechiba
1662
1663	ALLOCATE (veget_max(kjpindex,nvm),stat=ier)
1664	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget_max','','')
1665
1666	ALLOCATE (tot_bare_soil(kjpindex),stat=ier)
1667	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_bare_soil','','')
1668
1669	ALLOCATE (lai(kjpindex,nvm),stat=ier)
1670	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lai','','')
1671	lai(:,:)=undef_sechiba
1672
1673	ALLOCATE (frac_age(kjpindex,nvm,nleafages),stat=ier)
1674	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_age','','')
1675	frac_age(:,:,:)=undef_sechiba
1676
1677	ALLOCATE (height(kjpindex,nvm),stat=ier)
1678	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for height','','')
1679	height(:,:)=undef_sechiba
1680
1681	ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
1682	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_nobio','','')
1683	frac_nobio(:,:) = undef_sechiba
1684
1685	ALLOCATE (snow_nobio(kjpindex,nnobio),stat=ier)
1686	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio','','')
1687	snow_nobio(:,:) = undef_sechiba
1688
1689	ALLOCATE (snow_nobio_age(kjpindex,nnobio),stat=ier)
1690	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio_age','','')
1691	snow_nobio_age(:,:) = undef_sechiba
1692
1693	ALLOCATE (assim_param(kjpindex,nvm,npco2),stat=ier)
1694	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for assim_param','','')
1695
1696	!! 1.4 Initialize 1D array allocation
1697	ALLOCATE (vevapflo(kjpindex),stat=ier)
1698	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapflo','','')
1699	vevapflo(:)=zero
1700
1701	ALLOCATE (vevapsno(kjpindex),stat=ier)
1702	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapsno','','')
1703
1704	ALLOCATE (vevapnu(kjpindex),stat=ier)
1705	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapnu','','')
1706
1707	ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
1708	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for totfrac_nobio','','')
1709
1710	ALLOCATE (floodout(kjpindex),stat=ier)
1711	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for floodout','','')
1712
1713	ALLOCATE (runoff(kjpindex),stat=ier)
1714	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for runoff','','')
1715
1716	ALLOCATE (drainage(kjpindex),stat=ier)
1717	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drainage','','')
1718
1719	ALLOCATE (returnflow(kjpindex),stat=ier)
1720	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for returnflow','','')
1721	returnflow(:) = zero
1722
1723	ALLOCATE (reinfiltration(kjpindex),stat=ier)
1724	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinfiltration','','')
1725	reinfiltration(:) = zero
1726
1727	ALLOCATE (irrigation(kjpindex),stat=ier)
1728	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for irrigation','','')
1729	irrigation(:) = zero
1730
1731	ALLOCATE (z0h(kjpindex),stat=ier)
1732	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0h','','')
1733
1734	ALLOCATE (z0m(kjpindex),stat=ier)
1735	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0m','','')
1736
1737	ALLOCATE (roughheight(kjpindex),stat=ier)
1738	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for roughheight','','')
1739
1740	ALLOCATE (emis(kjpindex),stat=ier)
1741	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for emis','','')
1742
1743	ALLOCATE (tot_melt(kjpindex),stat=ier)
1744	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_melt','','')
1745
1746	ALLOCATE (vbeta(kjpindex),stat=ier)
1747	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta','','')
1748
1749	ALLOCATE (rau(kjpindex),stat=ier)
1750	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rau','','')
1751
1752	ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
1753	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for deadleaf_cover','','')
1754
1755	ALLOCATE (stempdiag(kjpindex, nslm),stat=ier)
1756	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for stempdiag','','')
1757
1758	ALLOCATE (ftempdiag(kjpindex, ngrnd),stat=ier)
1759	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ftempdiag','','')
1760
1761	ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
1762	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for co2_flux','','')
1763	co2_flux(:,:)=zero
1764
1765	ALLOCATE (shumdiag(kjpindex,nslm),stat=ier)
1766	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag','','')
1767
1768	ALLOCATE (shumdiag_perma(kjpindex,nslm),stat=ier)
1769	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag_perma','','')
1770
1771	ALLOCATE (litterhumdiag(kjpindex),stat=ier)
1772	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for litterhumdiag','','')
1773
1774	ALLOCATE (ptnlev1(kjpindex),stat=ier)
1775	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ptnlev1','','')
1776
1777	ALLOCATE (k_litt(kjpindex),stat=ier)
1778	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for k_litt','','')
1779
1780	!! 1.5 Initialize 2D array allocation
1781	ALLOCATE (vevapwet(kjpindex,nvm),stat=ier)
1782	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapwet','','')
1783	vevapwet(:,:)=undef_sechiba
1784
1785	ALLOCATE (transpir(kjpindex,nvm),stat=ier)
1786	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpir','','')
1787
1788	ALLOCATE (transpot(kjpindex,nvm),stat=ier)
1789	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpot','','')
1790
1791	ALLOCATE (qsintmax(kjpindex,nvm),stat=ier)
1792	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintmax','','')
1793
1794	ALLOCATE (rveget(kjpindex,nvm),stat=ier)
1795	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rveget','','')
1796
1797	ALLOCATE (rstruct(kjpindex,nvm),stat=ier)
1798	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rstruct','','')
1799
1800	ALLOCATE (pgflux(kjpindex),stat=ier)
1801	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for pgflux','','')
1802	pgflux(:)= 0.0
1803
1804	ALLOCATE (cgrnd_snow(kjpindex,nsnow),stat=ier)
1805	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cgrnd_snow','','')
1806	cgrnd_snow(:,:) = 0
1807
1808	ALLOCATE (dgrnd_snow(kjpindex,nsnow),stat=ier)
1809	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for dgrnd_snow','','')
1810	dgrnd_snow(:,:) = 0
1811
1812	ALLOCATE (lambda_snow(kjpindex),stat=ier)
1813	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lambda_snow','','')
1814	lambda_snow(:) = 0
1815
1816	ALLOCATE (temp_sol_add(kjpindex),stat=ier)
1817	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol_add','','')
1818
1819	ALLOCATE (gtemp(kjpindex),stat=ier)
1820	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gtemp','','')
1821
1822	ALLOCATE (frac_snow_veg(kjpindex),stat=ier)
1823	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_snow_veg','','')
1824
1825	ALLOCATE (frac_snow_nobio(kjpindex,nnobio),stat=ier)
1826	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_snow_nobio','','')
1827
1828	ALLOCATE (snowrho(kjpindex,nsnow),stat=ier)
1829	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowrho','','')
1830
1831	ALLOCATE (snowheat(kjpindex,nsnow),stat=ier)
1832	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowheat','','')
1833
1834	ALLOCATE (snowgrain(kjpindex,nsnow),stat=ier)
1835	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowgrain','','')
1836
1837	ALLOCATE (snowtemp(kjpindex,nsnow),stat=ier)
1838	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowtemp','','')
1839
1840	ALLOCATE (snowdz(kjpindex,nsnow),stat=ier)
1841	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowdz','','')
1842
1843	ALLOCATE (mc_layh(kjpindex, nslm),stat=ier)
1844	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mc_layh','','')
1845
1846	ALLOCATE (mcl_layh(kjpindex, nslm),stat=ier)
1847	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcl_layh','','')
1848
1849	ALLOCATE (soilmoist(kjpindex, nslm),stat=ier)
1850	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilmoist','','')
1851
1852
1853	!1.5 Irrigation related variables
1854	ALLOCATE (root_deficit(kjpindex),stat=ier)
1855	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for root_deficit','','') !
1856
1857	ALLOCATE (irrig_frac(kjpindex),stat=ier)
1858	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for irrig_frac','','') !
1859	irrigation(:) = zero
1860
1861	ALLOCATE (irrigated_next(kjpindex),stat=ier) !
1862	IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable irrigated_next','','') !
1863
1864	ALLOCATE (fraction_aeirrig_sw(kjpindex),stat=ier) !
1865	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fraction_aeirrig_sw','','')
1866
1867	!! 1.6 Initialize indexing table for the vegetation fields.
1868	! In SECHIBA we work on reduced grids but to store in the full 3D filed vegetation variable
1869	! we need another index table : indexveg, indexsoil, indexnobio and indexgrnd
1870	DO ji = 1, kjpindex
1871	!
1872	DO jv = 1, nlai+1
1873	indexlai((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1874	ENDDO
1875	!
1876	DO jv = 1, nvm
1877	indexveg((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1878	ENDDO
1879	!
1880	DO jv = 1, nstm
1881	indexsoil((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1882	ENDDO
1883	!
1884	DO jv = 1, nnobio
1885	indexnobio((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1886	ENDDO
1887	!
1888	DO jv = 1, ngrnd
1889	indexgrnd((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1890	ENDDO
1891	!
1892	DO jv = 1, nsnow
1893	indexsnow((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
1894	ENDDO
1895
1896	DO jv = 1, nslm
1897	indexnslm((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
1898	ENDDO
1899
1900	DO jv = 1, nslm
1901	indexlayer((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1902	ENDDO
1903	!
1904	DO jv = 1, 2
1905	indexalb((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1906	ENDDO
1907	!
1908	ENDDO
1909
1910	!! 2. Read the default value that will be put into variable which are not in the restart file
1911	CALL ioget_expval(val_exp)
1912
1913	IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_init done '
1914
1915	END SUBROUTINE sechiba_init
1916
1917
1918	!! ==============================================================================================================================\n
1919	!! SUBROUTINE : sechiba_clear
1920	!!
1921	!>\BRIEF Deallocate memory of sechiba's variables
1922	!!
1923	!! DESCRIPTION : None
1924	!!
1925	!! RECENT CHANGE(S): None
1926	!!
1927	!! MAIN OUTPUT VARIABLE(S): None
1928	!!
1929	!! REFERENCE(S) : None
1930	!!
1931	!! FLOWCHART : None
1932	!! \n
1933	!_ ================================================================================================================================
1934
1935	SUBROUTINE sechiba_clear()
1936
1937	!! 1. Initialize first run
1938
1939	l_first_sechiba=.TRUE.
1940
1941	!! 2. Deallocate dynamic variables of sechiba
1942
1943	IF ( ALLOCATED (indexveg)) DEALLOCATE (indexveg)
1944	IF ( ALLOCATED (indexlai)) DEALLOCATE (indexlai)
1945	IF ( ALLOCATED (indexsoil)) DEALLOCATE (indexsoil)
1946	IF ( ALLOCATED (indexnobio)) DEALLOCATE (indexnobio)
1947	IF ( ALLOCATED (indexsnow)) DEALLOCATE (indexsnow)
1948	IF ( ALLOCATED (indexgrnd)) DEALLOCATE (indexgrnd)
1949	IF ( ALLOCATED (indexlayer)) DEALLOCATE (indexlayer)
1950	IF ( ALLOCATED (indexnslm)) DEALLOCATE (indexnslm)
1951	IF ( ALLOCATED (indexalb)) DEALLOCATE (indexalb)
1952	IF ( ALLOCATED (flood_res)) DEALLOCATE (flood_res)
1953	IF ( ALLOCATED (flood_frac)) DEALLOCATE (flood_frac)
1954	IF ( ALLOCATED (snow)) DEALLOCATE (snow)
1955	IF ( ALLOCATED (snow_age)) DEALLOCATE (snow_age)
1956	IF ( ALLOCATED (drysoil_frac)) DEALLOCATE (drysoil_frac)
1957	IF ( ALLOCATED (evap_bare_lim)) DEALLOCATE (evap_bare_lim)
1958	IF ( ALLOCATED (evap_bare_lim_ns)) DEALLOCATE (evap_bare_lim_ns)
1959	IF ( ALLOCATED (evapot)) DEALLOCATE (evapot)
1960	IF ( ALLOCATED (evapot_corr)) DEALLOCATE (evapot_corr)
1961	IF ( ALLOCATED (humrel)) DEALLOCATE (humrel)
1962	IF ( ALLOCATED (vegstress)) DEALLOCATE (vegstress)
1963	IF ( ALLOCATED (soiltile)) DEALLOCATE (soiltile)
1964	IF ( ALLOCATED (fraclut)) DEALLOCATE (fraclut)
1965	IF ( ALLOCATED (nwdFraclut)) DEALLOCATE (nwdFraclut)
1966	IF ( ALLOCATED (njsc)) DEALLOCATE (njsc)
1967	IF ( ALLOCATED (reinf_slope)) DEALLOCATE (reinf_slope)
1968	IF ( ALLOCATED (reinf_slope_soil)) DEALLOCATE (reinf_slope_soil)
1969
1970	!salma: adding soil hydraulic params
1971	IF ( ALLOCATED (ks)) DEALLOCATE (ks)
1972	IF ( ALLOCATED (nvan)) DEALLOCATE (nvan)
1973	IF ( ALLOCATED (avan)) DEALLOCATE (avan)
1974	IF ( ALLOCATED (mcr)) DEALLOCATE (mcr)
1975	IF ( ALLOCATED (mcs)) DEALLOCATE (mcs)
1976	IF ( ALLOCATED (mcfc)) DEALLOCATE (mcfc)
1977	IF ( ALLOCATED (mcw)) DEALLOCATE (mcw)
1978	!end salma: adding soil hydraulic params
1979
1980	IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1)
1981	IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4)
1982	IF ( ALLOCATED (vbeta5)) DEALLOCATE (vbeta5)
1983	IF ( ALLOCATED (soilcap)) DEALLOCATE (soilcap)
1984	IF ( ALLOCATED (soilflx)) DEALLOCATE (soilflx)
1985	IF ( ALLOCATED (temp_sol)) DEALLOCATE (temp_sol)
1986	IF ( ALLOCATED (qsurf)) DEALLOCATE (qsurf)
1987	IF ( ALLOCATED (qsintveg)) DEALLOCATE (qsintveg)
1988	IF ( ALLOCATED (vbeta2)) DEALLOCATE (vbeta2)
1989	IF ( ALLOCATED (vbeta3)) DEALLOCATE (vbeta3)
1990	IF ( ALLOCATED (vbeta3pot)) DEALLOCATE (vbeta3pot)
1991	IF ( ALLOCATED (gsmean)) DEALLOCATE (gsmean)
1992	IF ( ALLOCATED (cimean)) DEALLOCATE (cimean)
1993	IF ( ALLOCATED (gpp)) DEALLOCATE (gpp)
1994	IF ( ALLOCATED (temp_growth)) DEALLOCATE (temp_growth)
1995	IF ( ALLOCATED (veget)) DEALLOCATE (veget)
1996	IF ( ALLOCATED (veget_max)) DEALLOCATE (veget_max)
1997	IF ( ALLOCATED (tot_bare_soil)) DEALLOCATE (tot_bare_soil)
1998	IF ( ALLOCATED (lai)) DEALLOCATE (lai)
1999	IF ( ALLOCATED (frac_age)) DEALLOCATE (frac_age)
2000	IF ( ALLOCATED (height)) DEALLOCATE (height)
2001	IF ( ALLOCATED (roughheight)) DEALLOCATE (roughheight)
2002	IF ( ALLOCATED (frac_nobio)) DEALLOCATE (frac_nobio)
2003	IF ( ALLOCATED (snow_nobio)) DEALLOCATE (snow_nobio)
2004	IF ( ALLOCATED (snow_nobio_age)) DEALLOCATE (snow_nobio_age)
2005	IF ( ALLOCATED (assim_param)) DEALLOCATE (assim_param)
2006	IF ( ALLOCATED (vevapflo)) DEALLOCATE (vevapflo)
2007	IF ( ALLOCATED (vevapsno)) DEALLOCATE (vevapsno)
2008	IF ( ALLOCATED (vevapnu)) DEALLOCATE (vevapnu)
2009	IF ( ALLOCATED (totfrac_nobio)) DEALLOCATE (totfrac_nobio)
2010	IF ( ALLOCATED (floodout)) DEALLOCATE (floodout)
2011	IF ( ALLOCATED (runoff)) DEALLOCATE (runoff)
2012	IF ( ALLOCATED (drainage)) DEALLOCATE (drainage)
2013	IF ( ALLOCATED (reinfiltration)) DEALLOCATE (reinfiltration)
2014	IF ( ALLOCATED (irrigation)) DEALLOCATE (irrigation)
2015	IF ( ALLOCATED (tot_melt)) DEALLOCATE (tot_melt)
2016	IF ( ALLOCATED (vbeta)) DEALLOCATE (vbeta)
2017	IF ( ALLOCATED (rau)) DEALLOCATE (rau)
2018	IF ( ALLOCATED (deadleaf_cover)) DEALLOCATE (deadleaf_cover)
2019	IF ( ALLOCATED (stempdiag)) DEALLOCATE (stempdiag)
2020	IF ( ALLOCATED (ftempdiag)) DEALLOCATE (ftempdiag)
2021	IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
2022	IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
2023	IF ( ALLOCATED (shumdiag_perma)) DEALLOCATE (shumdiag_perma)
2024	IF ( ALLOCATED (litterhumdiag)) DEALLOCATE (litterhumdiag)
2025	IF ( ALLOCATED (ptnlev1)) DEALLOCATE (ptnlev1)
2026	IF ( ALLOCATED (k_litt)) DEALLOCATE (k_litt)
2027	IF ( ALLOCATED (vevapwet)) DEALLOCATE (vevapwet)
2028	IF ( ALLOCATED (transpir)) DEALLOCATE (transpir)
2029	IF ( ALLOCATED (transpot)) DEALLOCATE (transpot)
2030	IF ( ALLOCATED (qsintmax)) DEALLOCATE (qsintmax)
2031	IF ( ALLOCATED (rveget)) DEALLOCATE (rveget)
2032	IF ( ALLOCATED (rstruct)) DEALLOCATE (rstruct)
2033	IF ( ALLOCATED (frac_snow_veg)) DEALLOCATE (frac_snow_veg)
2034	IF ( ALLOCATED (frac_snow_nobio)) DEALLOCATE (frac_snow_nobio)
2035	IF ( ALLOCATED (snowrho)) DEALLOCATE (snowrho)
2036	IF ( ALLOCATED (snowgrain)) DEALLOCATE (snowgrain)
2037	IF ( ALLOCATED (snowtemp)) DEALLOCATE (snowtemp)
2038	IF ( ALLOCATED (snowdz)) DEALLOCATE (snowdz)
2039	IF ( ALLOCATED (snowheat)) DEALLOCATE (snowheat)
2040	IF ( ALLOCATED (cgrnd_snow)) DEALLOCATE (cgrnd_snow)
2041	IF ( ALLOCATED (dgrnd_snow)) DEALLOCATE (dgrnd_snow)
2042	IF ( ALLOCATED (lambda_snow)) DEALLOCATE(lambda_snow)
2043	IF ( ALLOCATED (gtemp)) DEALLOCATE (gtemp)
2044	IF ( ALLOCATED (pgflux)) DEALLOCATE (pgflux)
2045	IF ( ALLOCATED (mc_layh)) DEALLOCATE (mc_layh)
2046	IF ( ALLOCATED (mcl_layh)) DEALLOCATE (mcl_layh)
2047	IF ( ALLOCATED (soilmoist)) DEALLOCATE (soilmoist)
2048	IF ( ALLOCATED (root_deficit)) DEALLOCATE (root_deficit)
2049	IF ( ALLOCATED (irrig_frac)) DEALLOCATE (irrig_frac)
2050	IF ( ALLOCATED (irrigated_next)) DEALLOCATE (irrigated_next)
2051
2052	!! 3. Clear all allocated memory
2053
2054	CALL pft_parameters_clear
2055	CALL slowproc_clear
2056	CALL diffuco_clear
2057	CALL enerbil_clear
2058	CALL hydrol_clear
2059	CALL thermosoil_clear
2060	CALL condveg_clear
2061	CALL routing_wrapper_clear
2062
2063	END SUBROUTINE sechiba_clear
2064
2065
2066	!! ==============================================================================================================================\n
2067	!! SUBROUTINE : sechiba_var_init
2068	!!
2069	!>\BRIEF Calculate air density as a function of air temperature and
2070	!! pressure for each terrestrial pixel.
2071	!!
2072	!! RECENT CHANGE(S): None
2073	!!
2074	!! MAIN OUTPUT VARIABLE(S): air density (::rau, kg m^{-3}).
2075	!!
2076	!! REFERENCE(S) : None
2077	!!
2078	!! FLOWCHART : None
2079	!! \n
2080	!_ ================================================================================================================================
2081
2082	SUBROUTINE sechiba_var_init (kjpindex, rau, pb, temp_air)
2083
2084	!! 0.1 Input variables
2085
2086	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
2087	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
2088	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
2089
2090	!! 0.2 Output variables
2091
2092	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rau !! Air density @tex $(kg m^{-3})$ @endtex
2093
2094	!! 0.3 Modified variables
2095
2096	!! 0.4 Local variables
2097
2098	INTEGER(i_std) :: ji !! Indices (unitless)
2099	!_ ================================================================================================================================
2100
2101	!! 1. Calculate intial air density (::rau)
2102
2103	DO ji = 1,kjpindex
2104	rau(ji) = pa_par_hpa * pb(ji) / (cte_molr*temp_air(ji))
2105	END DO
2106
2107	IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_var_init done '
2108
2109	END SUBROUTINE sechiba_var_init
2110
2111
2112	!! ==============================================================================================================================\n
2113	!! SUBROUTINE : sechiba_end
2114	!!
2115	!>\BRIEF Swap old for newly calculated soil temperature.
2116	!!
2117	!! RECENT CHANGE(S): None
2118	!!
2119	!! MAIN OUTPUT VARIABLE(S): soil temperature (::temp_sol; K)
2120	!!
2121	!! REFERENCE(S) : None
2122	!!
2123	!! FLOWCHART : None
2124	!! \n
2125	!! ================================================================================================================================
2126
2127	SUBROUTINE sechiba_end (kjpindex, temp_sol_new, temp_sol)
2128
2129
2130	!! 0.1 Input variables
2131
2132	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
2133	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature (K)
2134
2135	!! 0.2 Output variables
2136	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol !! Soil temperature (K)
2137
2138	!_ ================================================================================================================================
2139
2140	!! 1. Swap temperature
2141
2142	temp_sol(:) = temp_sol_new(:)
2143
2144	IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_end done '
2145
2146	END SUBROUTINE sechiba_end
2147
2148	!! ==============================================================================================================================\n
2149	!! SUBROUTINE : sechiba_interface_orchidee_inca
2150	!!
2151	!>\BRIEF make the interface between surface and atmospheric chemistry
2152	!!
2153	!! DESCRIPTION : This subroutine is called from INCA, the atmospheric chemistry model. It is used to transfer variables from ORCHIDEE to INCA.
2154	!!
2155	!! RECENT CHANGE(S): move from chemistry module to be more generic (feb - 2017)
2156	!!
2157	!! MAIN OUTPUT VARIABLE(S): emission COV to be transport by orchidee to inca in fields_out array
2158	!!
2159	!! REFERENCE(S) : None
2160	!!
2161	!! FLOWCHART : None
2162	!! \n
2163	!! ================================================================================================================================
2164	SUBROUTINE sechiba_interface_orchidee_inca( &
2165	nvm_out, veget_max_out, veget_frac_out, lai_out, snow_out, &
2166	field_out_COV_names, fields_out_COV, field_in_COV_names, fields_in_COV, &
2167	field_out_Nsoil_names, fields_out_Nsoil, nnspec_out)
2168
2169
2170	INTEGER, INTENT(out) :: nvm_out !! Number of vegetation types
2171	REAL(r_std), DIMENSION (:,:), INTENT(out) :: veget_max_out !! Max. fraction of vegetation type (LAI -> infty)
2172	REAL(r_std), DIMENSION (:,:), INTENT(out) :: veget_frac_out !! Fraction of vegetation type (unitless, 0-1)
2173	REAL(r_std), DIMENSION (:,:), INTENT(out) :: lai_out !! Surface foliere
2174	REAL(r_std), DIMENSION (:) , INTENT(out) :: snow_out !! Snow mass [Kg/m^2]
2175
2176	!
2177	! Optional arguments
2178	!
2179	! Names and fields for emission variables : to be transport by Orchidee to Inca
2180	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_out_COV_names
2181	REAL(r_std),DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: fields_out_COV
2182	!
2183	! Names and fields for deposit variables : to be transport from chemistry model by INCA to ORCHIDEE.
2184	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_in_COV_names
2185	REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(IN) :: fields_in_COV
2186
2187	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_out_Nsoil_names !! Not used before ORCHIDEE_3
2188	REAL(r_std),DIMENSION(:,:,:), OPTIONAL, INTENT(out) :: fields_out_Nsoil !! Not used before ORCHIDEE_3
2189	INTEGER, OPTIONAL, INTENT(out) :: nnspec_out !! Not used before ORCHIDEE_3
2190
2191
2192	IF (PRESENT(field_out_Nsoil_names) .OR. PRESENT(fields_out_Nsoil) .OR. PRESENT(nnspec_out)) THEN
2193	CALL ipslerr_p(3,'sechiba_interface_orchidee_inca','This version of Orchidee is not usable for coupling nitrogen with atmosphere',&
2194	'Please use Orchidee_3 or modify coupling between atm and surf','')
2195	ENDIF
2196
2197	! Variables always transmitted from sechiba to inca
2198	nvm_out = nvm
2199	veget_max_out(:,:) = veget_max(:,:)
2200	veget_frac_out(:,:) = veget(:,:)
2201	lai_out(:,:) = lai(:,:)
2202	snow_out(:) = snow(:)
2203
2204	! Call chemistry_flux_interface if at least one of variables field_out_names or
2205	! field_in_names is present in the argument list of sechiba_interface_orchidee_inca when called from inca.
2206	IF (PRESENT(field_out_COV_names) .AND. .NOT. PRESENT(field_in_COV_names)) THEN
2207	CALL chemistry_flux_interface(field_out_names=field_out_COV_names, fields_out=fields_out_COV)
2208	ELSE IF (.NOT. PRESENT(field_out_COV_names) .AND. PRESENT(field_in_COV_names)) THEN
2209	CALL chemistry_flux_interface(field_in_names=field_in_COV_names, fields_in=fields_in_COV)
2210	ELSE IF (PRESENT(field_out_COV_names) .AND. PRESENT(field_in_COV_names)) THEN
2211	CALL chemistry_flux_interface(field_out_names=field_out_COV_names, fields_out=fields_out_COV, &
2212	field_in_names=field_in_COV_names, fields_in=fields_in_COV)
2213	ENDIF
2214
2215	END SUBROUTINE sechiba_interface_orchidee_inca
2216
2217
2218	END MODULE sechiba
2219

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

Download in other formats: