Context Navigation

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

Last change on this file since 7796 was 7710, checked in by josefine.ghattas, 2 years ago
Integration of temperature of water in highres routing scheme. Done by Jan Polcer and Lucia Rinchiuso
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 116.9 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
521
522	!! 1.12 Initialize river routing
523	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
524
525	!! 1.12.1 Initialize river routing
526	CALL routing_wrapper_initialize( &
527	kjit, kjpindex, index, &
528	rest_id, hist_id, hist2_id, lalo, &
529	neighbours, resolution, contfrac, stempdiag, ftempdiag, &
530	soiltile, irrig_frac, veget_max, irrigated_next, &
531	returnflow, reinfiltration, irrigation, riverflow, &
532	coastalflow, flood_frac, flood_res )
533	ELSE
534	!! 1.12.2 No routing, set variables to zero
535	riverflow(:) = zero
536	coastalflow(:) = zero
537	returnflow(:) = zero
538	reinfiltration(:) = zero
539	irrigation(:) = zero
540	flood_frac(:) = zero
541	flood_res(:) = zero
542	ENDIF
543
544	!! 1.13 Write internal variables to output fields
545	z0m_out(:) = z0m(:)
546	z0h_out(:) = z0h(:)
547	emis_out(:) = emis(:)
548	qsurf_out(:) = qsurf(:)
549
550	!! 2. Output variables only once
551	zmaxh_glo(:) = zmaxh
552	CALL xios_orchidee_send_field("zmaxh",zmaxh_glo)
553
554	IF (printlev_loc>=3) WRITE(numout,*) 'sechiba_initialize done'
555
556	END SUBROUTINE sechiba_initialize
557
558	!! ==============================================================================================================================\n
559	!! SUBROUTINE : sechiba_main
560	!!
561	!>\BRIEF Main routine for the sechiba module performing three functions:
562	!! calculating temporal evolution of all variables and preparation of output and
563	!! restart files (during the last call only)
564	!!
565	!!\n DESCRIPTION : Main routine for the sechiba module.
566	!! One time step evolution consists of:
567	!! - call sechiba_var_init to do some initialization,
568	!! - call slowproc_main to do some daily calculations
569	!! - call diffuco_main for diffusion coefficient calculation,
570	!! - call enerbil_main for energy budget calculation,
571	!! - call hydrol_main for hydrologic processes calculation,
572	!! - call condveg_main for surface conditions such as roughness, albedo, and emmisivity,
573	!! - call thermosoil_main for soil thermodynamic calculation,
574	!! - call sechiba_end to swap previous to new fields.
575	!!
576	!! RECENT CHANGE(S): None
577	!!
578	!! MAIN OUTPUT VARIABLE(S): Hydrological variables (:: coastalflow and :: riverflow),
579	!! components of the energy budget (:: tsol_rad, :: vevapp, :: fluxsens,
580	!! :: temp_sol_new and :: fluxlat), surface characteristics (:: z0_out, :: emis_out,
581	!! :: tq_cdrag and :: albedo) and land use related CO2 fluxes (:: netco2flux and
582	!! :: fco2_lu, :: fco2_wh, ::fco2_ha)
583	!!
584	!! REFERENCE(S) :
585	!!
586	!! FLOWCHART :
587	!! \latexonly
588	!! \includegraphics[scale = 0.5]{sechibamainflow.png}
589	!! \endlatexonly
590	!! \n
591	!_ ================================================================================================================================
592
593	SUBROUTINE sechiba_main (kjit, kjpij, kjpindex, index, &
594	& ldrestart_read, ldrestart_write, &
595	& lalo, contfrac, neighbours, resolution,&
596	& zlev, u, v, qair, temp_air, epot_air, ccanopy, &
597	& tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
598	& precip_rain, precip_snow, lwdown, swnet, swdown, coszang, pb, &
599	& vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
600	& netco2flux, fco2_lu, fco2_wh, fco2_ha, &
601	& tsol_rad, temp_sol_new, qsurf_out, albedo, emis_out, z0m_out, z0h_out,&
602	& veget_out, lai_out, height_out, &
603	& rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
604
605	!! 0.1 Input variables
606
607	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
608	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
609	!! (unitless)
610	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
611	!! (unitless)
612	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
613	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (unitless)
614	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier (unitless)
615	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
616	!! (unitless)
617	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _History_ file identifier
618	!! (unitless)
619	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file
620	!! identifier (unitless)
621	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
622	!! (true/false)
623	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
624	!! (true/false)
625	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographic coordinates (latitude,longitude)
626	!! for grid cells (degrees)
627	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid
628	!! (unitless, 0-1)
629	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the pixels on the map.
630	!! Sechiba uses a reduced grid excluding oceans
631	!! ::index contains the indices of the
632	!! terrestrial pixels only! (unitless)
633	INTEGER(i_std), DIMENSION(kjpindex,NbNeighb), INTENT(in) :: neighbours !! Neighboring grid points if land!(unitless)
634	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size in x and y of the grid (m)
635
636	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed in direction u
637	!! @tex $(m.s^{-1})$ @endtex
638	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed in direction v
639	!! @tex $(m.s^{-1})$ @endtex
640	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
641	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
642	!! @tex $(kg kg^{-1})$ @endtex
643	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
644	!! @tex $(kg m^{-2})$ @endtex
645	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
646	!! @tex $(kg m^{-2})$ @endtex
647	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux
648	!! @tex $(W m^{-2})$ @endtex
649	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: coszang !! Cosine of the solar zenith angle (unitless)
650	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux
651	!! @tex $(W m^{-2})$ @endtex
652	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux
653	!! @tex $(W m^{-2})$ @endtex
654	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
655	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (??J)
656	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration in the canopy (ppm)
657	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! Coefficients A for T from the Planetary
658	!! Boundary Layer
659	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! Coefficients A for q from the Planetary
660	!! Boundary Layer
661	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! Coefficients B for T from the Planetary
662	!! Boundary Layer
663	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! Coefficients B for q from the Planetary
664	!! Boundary Layer
665	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
666
667
668	!! 0.2 Output variables
669
670	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: coastalflow !! Outflow on coastal points by small basins.
671	!! This is the water which flows in a disperse
672	!! way into the ocean
673	!! @tex $(kg dt_routing^{-1})$ @endtex
674	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: riverflow !! Outflow of the major rivers.
675	!! The flux will be located on the continental
676	!! grid but this should be a coastal point
677	!! @tex $(kg dt_routing^{-1})$ @endtex
678	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature
679	!! @tex $(W m^{-2})$ @endtex
680	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation
681	!! @tex $(kg m^{-2} days^{-1})$ @endtex
682
683	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf_out !! Surface specific humidity
684	!! @tex $(kg kg^{-1})$ @endtex
685	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m_out !! Surface roughness momentum (output diagnostic, m)
686	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h_out !! Surface roughness heat (output diagnostic, m)
687	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Surface albedo for visible and near-infrared
688	!! (unitless, 0-1)
689	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux
690	!! @tex $(W m^{-2})$ @endtex
691	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux
692	!! @tex $(W m^{-2})$ @endtex
693	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis_out !! Emissivity (unitless)
694	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: netco2flux !! Sum CO2 flux over PFTs
695	!! (gC/m2/dt_sechiba)
696	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_lu !! Land Cover Change CO2 flux (gC/m2/one_day)
697	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_wh !! Wood harvest CO2 flux (gC/m2/one_day)
698	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_ha !! Crop harvest CO2 flux (gC/m2/one_day)
699	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget_out !! Fraction of vegetation type (unitless, 0-1)
700	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: lai_out !! Leaf area index (m^2 m^{-2})
701	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height_out !! Vegetation Height (m)
702
703	!! 0.3 Modified
704
705	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: tq_cdrag !! Surface drag coefficient (-)
706	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! New ground temperature (K)
707
708	!! 0.4 local variables
709
710	INTEGER(i_std) :: ji, jv !! Index (unitless)
711	REAL(r_std), DIMENSION(kjpindex) :: histvar !! Computations for history files (unitless)
712	REAL(r_std), DIMENSION(kjpindex,nlut) :: histvar2 !! Computations for history files (unitless)
713	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
714	REAL(r_std), DIMENSION(kjpindex) :: sum_treefrac !! Total fraction occupied by trees (0-1, uniless)
715	REAL(r_std), DIMENSION(kjpindex) :: sum_grassfracC3 !! Total fraction occupied by C3 grasses (0-1, unitless)
716	REAL(r_std), DIMENSION(kjpindex) :: sum_grassfracC4 !! Total fraction occupied by C4 grasses (0-1, unitless)
717	REAL(r_std), DIMENSION(kjpindex) :: sum_cropfracC3 !! Total fraction occupied by C3 crops (0-1, unitess)
718	REAL(r_std), DIMENSION(kjpindex) :: sum_cropfracC4 !! Total fraction occupied by C4 crops (0-1, unitess)
719	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracNdlEvg!! Total fraction occupied by treeFracNdlEvg (0-1, unitess)
720	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracBdlEvg!! Total fraction occupied by treeFracBdlEvg (0-1, unitess)
721	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracNdlDcd!! Total fraction occupied by treeFracNdlDcd (0-1, unitess)
722	REAL(r_std), DIMENSION(kjpindex) :: sum_treeFracBdlDcd!! Total fraction occupied by treeFracBdlDcd (0-1, unitess)
723	REAL(r_std), DIMENSION(kjpindex) :: grndflux !! Net energy into soil (W/m2)
724	REAL(r_std), DIMENSION(kjpindex,nsnow) :: snowliq !! Liquid water content (m)
725	REAL(r_std), DIMENSION(kjpindex) :: snow_age_diag !! Only for diag, contains xios_default_val
726	REAL(r_std), DIMENSION(kjpindex,nnobio) :: snow_nobio_age_diag !! Only for diag, contains xios_default_val
727	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
728	REAL(r_std), DIMENSION(kjpindex,nlut) :: gpplut !! GPP on landuse tile, only for diagnostics
729
730
731	!_ ================================================================================================================================
732
733	IF (printlev_loc>=3) WRITE(numout,*) 'Start sechiba_main kjpindex =',kjpindex
734
735	!! 1. Initialize variables at each time step
736	CALL sechiba_var_init (kjpindex, rau, pb, temp_air)
737
738	!! 2. Compute diffusion coefficients
739	CALL diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, &
740	& zlev, z0m, z0h, roughheight, temp_sol, temp_air, temp_growth, rau, tq_cdrag, qsurf, qair, pb , &
741	& evap_bare_lim, evap_bare_lim_ns, evapot, evapot_corr, snow, flood_frac, flood_res, &
742	& frac_nobio, snow_nobio, totfrac_nobio, &
743	& swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
744	& vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, gsmean, rveget, rstruct, cimean, gpp, &
745	& lalo, neighbours, resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
746	& hist_id, hist2_id)
747
748	!! 3. Compute energy balance
749	CALL enerbil_main (kjit, kjpindex, &
750	& index, indexveg, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, &
751	& qair, peqAcoef, peqBcoef, pb, rau, vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
752	& emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
753	& vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo, temp_sol, tsol_rad, &
754	& temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id, &
755	& precip_rain, pgflux, snowdz, temp_sol_add)
756
757
758	!! 4. Compute hydrology
759	!! 4.1 Water balance from CWRR module (11 soil layers)
760	CALL hydrol_main (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjit, kjpindex, &
761	& index, indexveg, indexsoil, indexlayer, indexnslm, &
762	& temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, &
763	& qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age, &
764	& tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, &
765	& humrel, vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, evap_bare_lim_ns, flood_frac, flood_res, &
766	& shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope_soil,&
767	& rest_id, hist_id, hist2_id,&
768	& contfrac, stempdiag, &
769	& temp_air, pb, u, v, tq_cdrag, swnet, pgflux, &
770	& snowrho,snowtemp,snowgrain,snowdz,snowheat,snowliq, &
771	& grndflux,gtemp,tot_bare_soil, &
772	& lambda_snow,cgrnd_snow,dgrnd_snow,frac_snow_veg,temp_sol_add, &
773	& mc_layh, mcl_layh, soilmoist, root_deficit)
774
775
776	!! 6. Compute surface variables (emissivity, albedo and roughness)
777	CALL condveg_main (kjit, kjpindex, index, rest_id, hist_id, hist2_id, &
778	lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
779	zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
780	drysoil_frac, height, snowdz, snowrho, tot_bare_soil, &
781	temp_air, pb, u, v, lai, &
782	emis, albedo, z0m, z0h, roughheight, &
783	frac_snow_veg, frac_snow_nobio)
784
785
786	!! 7. Compute soil thermodynamics
787	CALL thermosoil_main (kjit, kjpindex, &
788	index, indexgrnd, mcs, &
789	temp_sol_new, snow, soilcap, soilflx, &
790	shumdiag_perma, stempdiag, ftempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
791	snowdz,snowrho,snowtemp,gtemp,pb,&
792	mc_layh, mcl_layh, soilmoist, njsc,frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
793	lambda_snow, cgrnd_snow, dgrnd_snow)
794
795
796	!! 8. Compute river routing
797	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
798	!! 8.1 River routing
799	CALL routing_wrapper_main (kjit, kjpindex, index, &
800	& lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
801	& drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, stempdiag, &
802	& ftempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id, &
803	& soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw)
804	ELSE
805	!! 8.2 No routing, set variables to zero
806	riverflow(:) = zero
807	coastalflow(:) = zero
808	returnflow(:) = zero
809	reinfiltration(:) = zero
810	irrigation(:) = zero
811	flood_frac(:) = zero
812	flood_res(:) = zero
813
814	CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba)
815	CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba)
816	ENDIF
817
818	!! 9. Compute slow processes (i.e. 'daily' and annual time step)
819	CALL slowproc_main (kjit, kjpij, kjpindex, &
820	index, indexveg, lalo, neighbours, resolution, contfrac, soiltile, fraclut, nwdFraclut, &
821	temp_air, temp_sol, stempdiag, &
822	vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
823	deadleaf_cover, &
824	assim_param, &
825	lai, frac_age, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
826	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
827	co2_flux, fco2_lu, fco2_wh, fco2_ha, temp_growth, tot_bare_soil, &
828	irrigated_next, irrig_frac, reinf_slope, reinf_slope_soil)
829
830
831	!! 9.2 Compute global CO2 flux
832	netco2flux(:) = zero
833	DO jv = 2,nvm
834	netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*(1-totfrac_nobio)
835	ENDDO
836
837	!! 10. Update the temperature (temp_sol) with newly computed values
838	CALL sechiba_end (kjpindex, temp_sol_new, temp_sol)
839
840
841	!! 11. Write internal variables to output fields
842	z0m_out(:) = z0m(:)
843	z0h_out(:) = z0h(:)
844	emis_out(:) = emis(:)
845	qsurf_out(:) = qsurf(:)
846	veget_out(:,:) = veget(:,:)
847	lai_out(:,:) = lai(:,:)
848	height_out(:,:) = height(:,:)
849
850	!! 12. Write global variables to history files
851	sum_treefrac(:) = zero
852	sum_grassfracC3(:) = zero
853	sum_grassfracC4(:) = zero
854	sum_cropfracC3(:) = zero
855	sum_cropfracC4(:) = zero
856	sum_treeFracNdlEvg(:) = zero
857	sum_treeFracBdlEvg(:) = zero
858	sum_treeFracNdlDcd(:) = zero
859	sum_treeFracBdlDcd(:) = zero
860
861	DO jv = 2, nvm
862	IF (is_tree(jv) .AND. natural(jv)) THEN
863	sum_treefrac(:) = sum_treefrac(:) + veget_max(:,jv)
864	ELSE IF ((.NOT. is_tree(jv)) .AND. natural(jv)) THEN
865	! Grass
866	IF (is_c4(jv)) THEN
867	sum_grassfracC4(:) = sum_grassfracC4(:) + veget_max(:,jv)
868	ELSE
869	sum_grassfracC3(:) = sum_grassfracC3(:) + veget_max(:,jv)
870	END IF
871	ELSE
872	! Crop and trees not natural
873	IF (is_c4(jv)) THEN
874	sum_cropfracC4(:) = sum_cropfracC4(:) + veget_max(:,jv)
875	ELSE
876	sum_cropfracC3(:) = sum_cropfracC3(:) + veget_max(:,jv)
877	END IF
878	ENDIF
879
880	IF (is_tree(jv)) THEN
881	IF (is_evergreen(jv)) THEN
882	IF (is_needleleaf(jv)) THEN
883	! Fraction for needleleaf evergreen trees (treeFracNdlEvg)
884	sum_treeFracNdlEvg(:) = sum_treeFracNdlEvg(:) + veget_max(:,jv)
885	ELSE
886	! Fraction for broadleaf evergreen trees (treeFracBdlEvg)
887	sum_treeFracBdlEvg(:) = sum_treeFracBdlEvg(:) + veget_max(:,jv)
888	END IF
889	ELSE IF (is_deciduous(jv)) THEN
890	IF (is_needleleaf(jv)) THEN
891	! Fraction for needleleaf deciduous trees (treeFracNdlDcd)
892	sum_treeFracNdlDcd(:) = sum_treeFracNdlDcd(:) + veget_max(:,jv)
893	ELSE
894	! Fraction for broadleafs deciduous trees (treeFracBdlDcd)
895	sum_treeFracBdlDcd(:) = sum_treeFracBdlDcd(:) + veget_max(:,jv)
896	END IF
897	END IF
898	END IF
899	ENDDO
900
901	histvar(:)=zero
902	DO jv = 2, nvm
903	IF (is_deciduous(jv)) THEN
904	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
905	ENDIF
906	ENDDO
907	CALL xios_orchidee_send_field("treeFracPrimDec",histvar)
908
909	histvar(:)=zero
910	DO jv = 2, nvm
911	IF (is_evergreen(jv)) THEN
912	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
913	ENDIF
914	ENDDO
915	CALL xios_orchidee_send_field("treeFracPrimEver",histvar)
916
917	histvar(:)=zero
918	DO jv = 2, nvm
919	IF ( .NOT.(is_c4(jv)) ) THEN
920	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
921	ENDIF
922	ENDDO
923	CALL xios_orchidee_send_field("c3PftFrac",histvar)
924
925	histvar(:)=zero
926	DO jv = 2, nvm
927	IF ( is_c4(jv) ) THEN
928	histvar(:) = histvar(:) + veget_max(:,jv)100contfrac
929	ENDIF
930	ENDDO
931	CALL xios_orchidee_send_field("c4PftFrac",histvar)
932
933	CALL xios_orchidee_send_field("temp_sol_new",temp_sol_new)
934	CALL xios_orchidee_send_field("fluxsens",fluxsens)
935	CALL xios_orchidee_send_field("fluxlat",fluxlat)
936
937
938	! Add XIOS default value where no snow
939	DO ji=1,kjpindex
940	IF (snow(ji) .GT. zero) THEN
941	snow_age_diag(ji) = snow_age(ji)
942	snow_nobio_age_diag(ji,:) = snow_nobio_age(ji,:)
943
944	snowage_glob(ji) = snow_age(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji)) + &
945	SUM(snow_nobio_age(ji,:)frac_snow_nobio(ji,:)frac_nobio(ji,:))
946	IF (snowage_glob(ji) .NE. 0) snowage_glob(ji) = snowage_glob(ji) / &
947	(frac_snow_veg(ji)(1-totfrac_nobio(ji)) + SUM(frac_snow_nobio(ji,:)frac_nobio(ji,:)))
948	ELSE
949	snow_age_diag(ji) = xios_default_val
950	snow_nobio_age_diag(ji,:) = xios_default_val
951	snowage_glob(ji) = xios_default_val
952	END IF
953	END DO
954
955	CALL xios_orchidee_send_field("snow",snow)
956	CALL xios_orchidee_send_field("snowage",snow_age_diag)
957	CALL xios_orchidee_send_field("snownobio",snow_nobio)
958	CALL xios_orchidee_send_field("snownobioage",snow_nobio_age_diag)
959	CALL xios_orchidee_send_field("snowage_glob",snowage_glob)
960
961	CALL xios_orchidee_send_field("frac_snow", SUM(frac_snow_nobio,2)totfrac_nobio+frac_snow_veg(1-totfrac_nobio))
962	CALL xios_orchidee_send_field("frac_snow_veg", frac_snow_veg)
963	CALL xios_orchidee_send_field("frac_snow_nobio", frac_snow_nobio)
964	CALL xios_orchidee_send_field("reinf_slope",reinf_slope)
965	CALL xios_orchidee_send_field("njsc",REAL(njsc, r_std))
966	CALL xios_orchidee_send_field("vegetfrac",veget)
967	CALL xios_orchidee_send_field("maxvegetfrac",veget_max)
968	CALL xios_orchidee_send_field("irrigmap_dyn",irrigated_next)
969	CALL xios_orchidee_send_field("aei_sw",fraction_aeirrig_sw)
970	CALL xios_orchidee_send_field("nobiofrac",frac_nobio)
971	CALL xios_orchidee_send_field("soiltile",soiltile)
972	CALL xios_orchidee_send_field("rstruct",rstruct)
973	CALL xios_orchidee_send_field("gpp",gpp/dt_sechiba)
974	CALL xios_orchidee_send_field("gpp_ipcc2",SUM(gpp,dim=2)/dt_sechiba)
975
976	histvar(:)=zero
977	DO jv = 2, nvm
978	IF ( .NOT. is_tree(jv) .AND. natural(jv) ) THEN
979	histvar(:) = histvar(:) + gpp(:,jv)
980	ENDIF
981	ENDDO
982	CALL xios_orchidee_send_field("gppgrass",histvar/dt_sechiba)
983
984	histvar(:)=zero
985	DO jv = 2, nvm
986	IF ( (.NOT. is_tree(jv)) .AND. (.NOT. natural(jv)) ) THEN
987	histvar(:) = histvar(:) + gpp(:,jv)
988	ENDIF
989	ENDDO
990	CALL xios_orchidee_send_field("gppcrop",histvar/dt_sechiba)
991
992	histvar(:)=zero
993	DO jv = 2, nvm
994	IF ( is_tree(jv) ) THEN
995	histvar(:) = histvar(:) + gpp(:,jv)
996	ENDIF
997	ENDDO
998	CALL xios_orchidee_send_field("gpptree",histvar/dt_sechiba)
999	CALL xios_orchidee_send_field("nee",co2_flux/1.e3/one_day)
1000	CALL xios_orchidee_send_field("drysoil_frac",drysoil_frac)
1001	CALL xios_orchidee_send_field("vevapflo",vevapflo/dt_sechiba)
1002	CALL xios_orchidee_send_field("k_litt",k_litt)
1003	CALL xios_orchidee_send_field("beta",vbeta)
1004	CALL xios_orchidee_send_field("vbeta1",vbeta1)
1005	CALL xios_orchidee_send_field("vbeta2",vbeta2)
1006	CALL xios_orchidee_send_field("vbeta3",vbeta3)
1007	CALL xios_orchidee_send_field("vbeta4",vbeta4)
1008	CALL xios_orchidee_send_field("vbeta5",vbeta5)
1009	CALL xios_orchidee_send_field("gsmean",gsmean)
1010	CALL xios_orchidee_send_field("cimean",cimean)
1011	CALL xios_orchidee_send_field("rveget",rveget)
1012
1013	histvar(:)=SUM(vevapwet(:,:),dim=2)
1014	CALL xios_orchidee_send_field("evspsblveg",histvar/dt_sechiba)
1015	histvar(:)= vevapnu(:)+vevapsno(:)
1016	CALL xios_orchidee_send_field("evspsblsoi",histvar/dt_sechiba)
1017	histvar(:)=SUM(transpir(:,:),dim=2)
1018	CALL xios_orchidee_send_field("tran",histvar/dt_sechiba)
1019
1020	! For CMIP6 data request: the following fractions are fractions of the total grid-cell,
1021	! which explains the multiplication by contfrac
1022	CALL xios_orchidee_send_field("treeFrac",sum_treefrac(:)100contfrac(:))
1023	CALL xios_orchidee_send_field("grassFracC3",sum_grassfracC3(:)100contfrac(:))
1024	CALL xios_orchidee_send_field("grassFracC4",sum_grassfracC4(:)100contfrac(:))
1025	CALL xios_orchidee_send_field("cropFracC3",sum_cropfracC3(:)100contfrac(:))
1026	CALL xios_orchidee_send_field("cropFracC4",sum_cropfracC4(:)100contfrac(:))
1027	CALL xios_orchidee_send_field("treeFracNdlEvg",sum_treeFracNdlEvg(:)100contfrac(:))
1028	CALL xios_orchidee_send_field("treeFracBdlEvg",sum_treeFracBdlEvg(:)100contfrac(:))
1029	CALL xios_orchidee_send_field("treeFracNdlDcd",sum_treeFracNdlDcd(:)100contfrac(:))
1030	CALL xios_orchidee_send_field("treeFracBdlDcd",sum_treeFracBdlDcd(:)100contfrac(:))
1031
1032	histvar(:)=veget_max(:,1)100contfrac(:)
1033	CALL xios_orchidee_send_field("baresoilFrac",histvar)
1034	histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)100contfrac(:)
1035	CALL xios_orchidee_send_field("residualFrac",histvar)
1036
1037	! For CMIP6 data request: cnc = canopy cover fraction over land area
1038	histvar(:)=zero
1039	DO jv=2,nvm
1040	histvar(:) = histvar(:) + veget_max(:,jv)*100
1041	END DO
1042	CALL xios_orchidee_send_field("cnc",histvar)
1043
1044	CALL xios_orchidee_send_field("tsol_rad",tsol_rad-273.15)
1045	CALL xios_orchidee_send_field("qsurf",qsurf)
1046	CALL xios_orchidee_send_field("emis",emis)
1047	CALL xios_orchidee_send_field("z0m",z0m)
1048	CALL xios_orchidee_send_field("z0h",z0h)
1049	CALL xios_orchidee_send_field("roughheight",roughheight)
1050	CALL xios_orchidee_send_field("lai",lai)
1051	histvar(:)=zero
1052	DO ji = 1, kjpindex
1053	IF (SUM(veget_max(ji,:)) > zero) THEN
1054	DO jv=2,nvm
1055	histvar(ji) = histvar(ji) + veget_max(ji,jv)*lai(ji,jv)/SUM(veget_max(ji,:))
1056	END DO
1057	END IF
1058	END DO
1059	CALL xios_orchidee_send_field("LAImean",histvar)
1060	CALL xios_orchidee_send_field("vevapsno",vevapsno/dt_sechiba)
1061	CALL xios_orchidee_send_field("vevapp",vevapp/dt_sechiba)
1062	CALL xios_orchidee_send_field("vevapnu",vevapnu/dt_sechiba)
1063	CALL xios_orchidee_send_field("transpir",transpir*one_day/dt_sechiba)
1064	CALL xios_orchidee_send_field("transpot",transpot*one_day/dt_sechiba)
1065	CALL xios_orchidee_send_field("inter",vevapwet*one_day/dt_sechiba)
1066	histvar(:)=zero
1067	DO jv=1,nvm
1068	histvar(:) = histvar(:) + vevapwet(:,jv)
1069	ENDDO
1070	CALL xios_orchidee_send_field("ECanop",histvar/dt_sechiba)
1071	histvar(:)=zero
1072	DO jv=1,nvm
1073	histvar(:) = histvar(:) + transpir(:,jv)
1074	ENDDO
1075	CALL xios_orchidee_send_field("TVeg",histvar/dt_sechiba)
1076
1077
1078	!! Calculate diagnostic variables on Landuse tiles for LUMIP/CMIP6
1079
1080	! Calculate fraction of landuse tiles related to the whole grid cell
1081	DO jv=1,nlut
1082	histvar2(:,jv) = fraclut(:,jv) * contfrac(:)
1083	END DO
1084	CALL xios_orchidee_send_field("fraclut",histvar2)
1085
1086	CALL xios_orchidee_send_field("nwdFraclut",nwdFraclut(:,:))
1087
1088	! Calculate GPP on landuse tiles
1089	! val_exp is used as missing value where the values are not known i.e. where the tile is not represented
1090	! or for pasture (id_pst) or urban land (id_urb).
1091	gpplut(:,:)=0
1092	DO jv=1,nvm
1093	IF (natural(jv)) THEN
1094	gpplut(:,id_psl) = gpplut(:,id_psl) + gpp(:,jv)
1095	ELSE
1096	gpplut(:,id_crp) = gpplut(:,id_crp) + gpp(:,jv)
1097	ENDIF
1098	END DO
1099
1100	! Transform from gC/m2/s into kgC/m2/s
1101	WHERE (fraclut(:,id_psl)>min_sechiba)
1102	gpplut(:,id_psl) = gpplut(:,id_psl)/fraclut(:,id_psl)/1000
1103	ELSEWHERE
1104	gpplut(:,id_psl) = xios_default_val
1105	END WHERE
1106	WHERE (fraclut(:,id_crp)>min_sechiba)
1107	gpplut(:,id_crp) = gpplut(:,id_crp)/fraclut(:,id_crp)/1000
1108	ELSEWHERE
1109	gpplut(:,id_crp) = xios_default_val
1110	END WHERE
1111	gpplut(:,id_pst) = xios_default_val
1112	gpplut(:,id_urb) = xios_default_val
1113
1114	CALL xios_orchidee_send_field("gpplut",gpplut)
1115
1116
1117	IF ( .NOT. almaoutput ) THEN
1118	! Write history file in IPSL-format
1119	CALL histwrite_p(hist_id, 'beta', kjit, vbeta, kjpindex, index)
1120	CALL histwrite_p(hist_id, 'z0m', kjit, z0m, kjpindex, index)
1121	CALL histwrite_p(hist_id, 'z0h', kjit, z0h, kjpindex, index)
1122	CALL histwrite_p(hist_id, 'roughheight', kjit, roughheight, kjpindex, index)
1123	CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1124	CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1125	CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1126	CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1127	CALL histwrite_p(hist_id, 'subli', kjit, vevapsno, kjpindex, index)
1128	CALL histwrite_p(hist_id, 'evapnu', kjit, vevapnu, kjpindex, index)
1129	CALL histwrite_p(hist_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
1130	CALL histwrite_p(hist_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
1131	CALL histwrite_p(hist_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
1132	CALL histwrite_p(hist_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
1133	CALL histwrite_p(hist_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
1134	CALL histwrite_p(hist_id, 'vbeta4', kjit, vbeta4, kjpindex, index)
1135	CALL histwrite_p(hist_id, 'vbeta5', kjit, vbeta5, kjpindex, index)
1136	CALL histwrite_p(hist_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
1137	CALL histwrite_p(hist_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
1138	CALL histwrite_p(hist_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
1139	CALL histwrite_p(hist_id, 'snow', kjit, snow, kjpindex, index)
1140	CALL histwrite_p(hist_id, 'snowage', kjit, snow_age, kjpindex, index)
1141	CALL histwrite_p(hist_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
1142	CALL histwrite_p(hist_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
1143
1144	CALL histwrite_p(hist_id, 'grndflux', kjit, grndflux, kjpindex,index)
1145	CALL histwrite_p(hist_id, 'snowtemp',kjit,snowtemp,kjpindex*nsnow,indexsnow)
1146	CALL histwrite_p(hist_id, 'snowliq', kjit,snowliq,kjpindex*nsnow,indexsnow)
1147	CALL histwrite_p(hist_id, 'snowdz', kjit,snowdz,kjpindex*nsnow,indexsnow)
1148	CALL histwrite_p(hist_id, 'snowrho', kjit,snowrho,kjpindex*nsnow,indexsnow)
1149	CALL histwrite_p(hist_id, 'snowgrain',kjit,snowgrain,kjpindex*nsnow,indexsnow)
1150	CALL histwrite_p(hist_id, 'snowheat',kjit,snowheat,kjpindex*nsnow,indexsnow)
1151
1152	CALL histwrite_p(hist_id,'frac_snow_veg',kjit,frac_snow_veg,kjpindex,index)
1153	CALL histwrite_p(hist_id, 'frac_snow_nobio', kjit,frac_snow_nobio,kjpindex*nnobio, indexnobio)
1154	CALL histwrite_p(hist_id, 'pgflux',kjit,pgflux,kjpindex,index)
1155	CALL histwrite_p(hist_id, 'soiltile', kjit, soiltile, kjpindex*nstm, indexsoil)
1156
1157	CALL histwrite_p(hist_id, 'soilindex', kjit, REAL(njsc, r_std), kjpindex, index)
1158	CALL histwrite_p(hist_id, 'reinf_slope', kjit, reinf_slope, kjpindex, index)
1159	CALL histwrite_p(hist_id, 'k_litt', kjit, k_litt, kjpindex, index)
1160
1161	IF ( do_floodplains ) THEN
1162	CALL histwrite_p(hist_id, 'evapflo', kjit, vevapflo, kjpindex, index)
1163	CALL histwrite_p(hist_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
1164	ENDIF
1165
1166	CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)
1167	CALL histwrite_p(hist_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
1168	CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
1169
1170	IF ( ok_stomate ) THEN
1171	CALL histwrite_p(hist_id, 'nee', kjit, co2_flux/1.e3/one_day, kjpindex*nvm, indexveg)
1172	ENDIF
1173
1174	histvar(:)=SUM(vevapwet(:,:),dim=2)
1175	CALL histwrite_p(hist_id, 'evspsblveg', kjit, histvar, kjpindex, index)
1176
1177	histvar(:)= vevapnu(:)+vevapsno(:)
1178	CALL histwrite_p(hist_id, 'evspsblsoi', kjit, histvar, kjpindex, index)
1179
1180	histvar(:)=SUM(transpir(:,:),dim=2)
1181	CALL histwrite_p(hist_id, 'tran', kjit, histvar, kjpindex, index)
1182
1183	histvar(:)= sum_treefrac(:)100contfrac(:)
1184	CALL histwrite_p(hist_id, 'treeFrac', kjit, histvar, kjpindex, index)
1185
1186	histvar(:)= (sum_grassfracC3(:)+sum_grassfracC4(:))100contfrac(:)
1187	CALL histwrite_p(hist_id, 'grassFrac', kjit, histvar, kjpindex, index)
1188
1189	histvar(:)= (sum_cropfracC3(:)+sum_cropfracC4(:))100contfrac(:)
1190	CALL histwrite_p(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)
1191
1192	histvar(:)=veget_max(:,1)100contfrac(:)
1193	CALL histwrite_p(hist_id, 'baresoilFrac', kjit, histvar, kjpindex, index)
1194
1195	histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)100contfrac(:)
1196	CALL histwrite_p(hist_id, 'residualFrac', kjit, histvar, kjpindex, index)
1197	ELSE
1198	! Write history file in ALMA format
1199	CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1200	CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1201	CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1202	CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1203	CALL histwrite_p(hist_id, 'ESoil', kjit, vevapnu, kjpindex, index)
1204	CALL histwrite_p(hist_id, 'EWater', kjit, vevapflo, kjpindex, index)
1205	CALL histwrite_p(hist_id, 'SWE', kjit, snow, kjpindex, index)
1206	histvar(:)=zero
1207	DO jv=1,nvm
1208	histvar(:) = histvar(:) + transpir(:,jv)
1209	ENDDO
1210	CALL histwrite_p(hist_id, 'TVeg', kjit, histvar, kjpindex, index)
1211	histvar(:)=zero
1212	DO jv=1,nvm
1213	histvar(:) = histvar(:) + vevapwet(:,jv)
1214	ENDDO
1215	CALL histwrite_p(hist_id, 'ECanop', kjit, histvar, kjpindex, index)
1216	CALL histwrite_p(hist_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
1217	!
1218	CALL histwrite_p(hist_id, 'Z0m', kjit, z0m, kjpindex, index)
1219	CALL histwrite_p(hist_id, 'Z0h', kjit, z0h, kjpindex, index)
1220	CALL histwrite_p(hist_id, 'EffectHeight', kjit, roughheight, kjpindex, index)
1221	!
1222	IF ( do_floodplains ) THEN
1223	CALL histwrite_p(hist_id, 'Qflood', kjit, vevapflo, kjpindex, index)
1224	CALL histwrite_p(hist_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
1225	ENDIF
1226
1227	CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)
1228	CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
1229	CALL histwrite_p(hist_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
1230
1231	IF ( ok_stomate ) THEN
1232	CALL histwrite_p(hist_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)
1233	ENDIF
1234	ENDIF ! almaoutput
1235
1236	!! 13. Write additional output file with higher frequency
1237	IF ( hist2_id > 0 ) THEN
1238	IF ( .NOT. almaoutput ) THEN
1239	! Write history file in IPSL-format
1240	CALL histwrite_p(hist2_id, 'tsol_rad', kjit, tsol_rad, kjpindex, index)
1241	CALL histwrite_p(hist2_id, 'qsurf', kjit, qsurf, kjpindex, index)
1242	CALL histwrite_p(hist2_id, 'albedo', kjit, albedo, kjpindex*2, indexalb)
1243	CALL histwrite_p(hist2_id, 'emis', kjit, emis, kjpindex, index)
1244	CALL histwrite_p(hist2_id, 'beta', kjit, vbeta, kjpindex, index)
1245	CALL histwrite_p(hist2_id, 'z0m', kjit, z0m, kjpindex, index)
1246	CALL histwrite_p(hist2_id, 'z0h', kjit, z0h, kjpindex, index)
1247	CALL histwrite_p(hist2_id, 'roughheight', kjit, roughheight, kjpindex, index)
1248	CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1249	CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1250	CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1251	CALL histwrite_p(hist2_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1252	CALL histwrite_p(hist2_id, 'subli', kjit, vevapsno, kjpindex, index)
1253	IF ( do_floodplains ) THEN
1254	CALL histwrite_p(hist2_id, 'vevapflo', kjit, vevapflo, kjpindex, index)
1255	CALL histwrite_p(hist2_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
1256	ENDIF
1257	CALL histwrite_p(hist2_id, 'vevapnu', kjit, vevapnu, kjpindex, index)
1258	CALL histwrite_p(hist2_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
1259	CALL histwrite_p(hist2_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
1260	CALL histwrite_p(hist2_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
1261	CALL histwrite_p(hist2_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
1262	CALL histwrite_p(hist2_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
1263	CALL histwrite_p(hist2_id, 'vbeta4', kjit, vbeta4, kjpindex, index)
1264	CALL histwrite_p(hist2_id, 'vbeta5', kjit, vbeta5, kjpindex, index)
1265	CALL histwrite_p(hist2_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
1266	CALL histwrite_p(hist2_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
1267	CALL histwrite_p(hist2_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
1268	CALL histwrite_p(hist2_id, 'snow', kjit, snow, kjpindex, index)
1269	CALL histwrite_p(hist2_id, 'snowage', kjit, snow_age, kjpindex, index)
1270	CALL histwrite_p(hist2_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
1271	CALL histwrite_p(hist2_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
1272	CALL histwrite_p(hist2_id, 'soilindex', kjit, REAL(njsc, r_std), kjpindex, index)
1273	CALL histwrite_p(hist2_id, 'reinf_slope', kjit, reinf_slope, kjpindex, index)
1274
1275	CALL histwrite_p(hist2_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)
1276	CALL histwrite_p(hist2_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
1277	CALL histwrite_p(hist2_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
1278
1279	IF ( ok_stomate ) THEN
1280	CALL histwrite_p(hist2_id, 'nee', kjit, co2_flux/1.e3/one_day, kjpindex*nvm, indexveg)
1281	ENDIF
1282	ELSE
1283	! Write history file in ALMA format
1284	CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1285	CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1286	CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1287	CALL histwrite_p(hist2_id, 'ESoil', kjit, vevapnu, kjpindex, index)
1288	IF ( do_floodplains ) THEN
1289	CALL histwrite_p(hist2_id, 'EWater', kjit, vevapflo, kjpindex, index)
1290	CALL histwrite_p(hist2_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
1291	ENDIF
1292	CALL histwrite_p(hist2_id, 'SWE', kjit, snow, kjpindex, index)
1293	histvar(:)=zero
1294	DO jv=1,nvm
1295	histvar(:) = histvar(:) + transpir(:,jv)
1296	ENDDO
1297	CALL histwrite_p(hist2_id, 'TVeg', kjit, histvar, kjpindex, index)
1298	histvar(:)=zero
1299	DO jv=1,nvm
1300	histvar(:) = histvar(:) + vevapwet(:,jv)
1301	ENDDO
1302	CALL histwrite_p(hist2_id, 'ECanop', kjit, histvar, kjpindex, index)
1303	CALL histwrite_p(hist2_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
1304	CALL histwrite_p(hist2_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
1305
1306	IF ( ok_stomate ) THEN
1307	CALL histwrite_p(hist2_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)
1308	ENDIF
1309	ENDIF ! almaoutput
1310	ENDIF ! hist2_id
1311
1312
1313	!! Change the vegetation fractions if a new map was read in slowproc. This is done
1314	!! after lcchange has been done in stomatelpj
1315	IF (done_stomate_lcchange) THEN
1316	CALL slowproc_change_frac(kjpindex, lai, &
1317	veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile, fraclut, nwdFraclut)
1318	done_stomate_lcchange=.FALSE.
1319	END IF
1320
1321	!! 14. If it is the last time step, write restart files
1322	IF (ldrestart_write) THEN
1323	CALL sechiba_finalize( &
1324	kjit, kjpij, kjpindex, index, rest_id, &
1325	tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad)
1326	END IF
1327
1328	END SUBROUTINE sechiba_main
1329
1330
1331	!! =============================================================================================================================
1332	!! SUBROUTINE: sechiba_finalize
1333	!!
1334	!>\BRIEF Finalize all modules by calling their "_finalize" subroutines.
1335	!!
1336	!! DESCRIPTION: Finalize all modules by calling their "_finalize" subroutines. These subroutines will write variables to
1337	!! restart file.
1338	!!
1339	!! \n
1340	!_ ==============================================================================================================================
1341
1342	SUBROUTINE sechiba_finalize( &
1343	kjit, kjpij, kjpindex, index, rest_id, &
1344	tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad)
1345
1346	!! 0.1 Input variables
1347	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
1348	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
1349	!! (unitless)
1350	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
1351	!! (unitless)
1352	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
1353	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the pixels on the map.
1354	!! Sechiba uses a reduced grid excluding oceans
1355	!! ::index contains the indices of the
1356	!! terrestrial pixels only! (unitless)
1357	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tsol_rad !! Radiative surface temperature
1358	!! @tex $(W m^{-2})$ @endtex
1359	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapp !! Total of evaporation
1360	!! @tex $(kg m^{-2} days^{-1})$ @endtex
1361	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxsens !! Sensible heat flux
1362	!! @tex $(W m^{-2})$ @endtex
1363	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxlat !! Latent heat flux
1364	!! @tex $(W m^{-2})$ @endtex
1365	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tq_cdrag !! Surface drag coefficient (-)
1366
1367	!! 0.2 Local variables
1368	INTEGER(i_std) :: ji, jv !! Index (unitless)
1369	REAL(r_std), DIMENSION(kjpindex) :: histvar !! Computations for history files (unitless)
1370	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
1371
1372
1373	!! Write restart file for the next simulation from SECHIBA and other modules
1374
1375	IF (printlev_loc>=3) WRITE (numout,*) 'Start sechiba_finalize for writing restart files'
1376
1377	!! 1. Call diffuco_finalize to write restart files
1378	CALL diffuco_finalize (kjit, kjpindex, rest_id, rstruct )
1379
1380	!! 2. Call energy budget to write restart files
1381	CALL enerbil_finalize (kjit, kjpindex, rest_id, &
1382	evapot, evapot_corr, temp_sol, tsol_rad, &
1383	qsurf, fluxsens, fluxlat, vevapp )
1384
1385	!! 3. Call hydrology to write restart files
1386	CALL hydrol_finalize( kjit, kjpindex, rest_id, vegstress, &
1387	qsintveg, humrel, snow, snow_age, snow_nobio, &
1388	snow_nobio_age, snowrho, snowtemp, snowdz, &
1389	snowheat, snowgrain, &
1390	drysoil_frac, evap_bare_lim, evap_bare_lim_ns)
1391
1392	!! 4. Call condveg to write surface variables to restart files
1393	CALL condveg_finalize (kjit, kjpindex, rest_id, z0m, z0h, roughheight)
1394
1395	!! 5. Call soil thermodynamic to write restart files
1396	CALL thermosoil_finalize (kjit, kjpindex, rest_id, gtemp, &
1397	soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
1398
1399
1400	!! 6. Add river routing to restart files
1401	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
1402	!! 6.1 Call river routing to write restart files
1403	CALL routing_wrapper_finalize( kjit, kjpindex, rest_id, flood_frac, flood_res )
1404	ELSE
1405	!! 6.2 No routing, set variables to zero
1406	reinfiltration(:) = zero
1407	returnflow(:) = zero
1408	irrigation(:) = zero
1409	flood_frac(:) = zero
1410	flood_res(:) = zero
1411	ENDIF
1412
1413	!! 7. Call slowproc_main to add 'daily' and annual variables to restart file
1414	CALL slowproc_finalize (kjit, kjpindex, rest_id, index, &
1415	njsc, lai, height, veget, &
1416	frac_nobio, veget_max, reinf_slope, &
1417	ks, nvan, avan, mcr, &
1418	mcs, mcfc, mcw, &
1419	assim_param, frac_age, fraction_aeirrig_sw)
1420
1421	IF (printlev_loc>=3) WRITE (numout,*) 'sechiba_finalize done'
1422
1423	END SUBROUTINE sechiba_finalize
1424
1425
1426	!! ==============================================================================================================================\n
1427	!! SUBROUTINE : sechiba_init
1428	!!
1429	!>\BRIEF Dynamic allocation of the variables, the dimensions of the
1430	!! variables are determined by user-specified settings.
1431	!!
1432	!! DESCRIPTION : The domain size (:: kjpindex) is used to allocate the correct
1433	!! dimensions to all variables in sechiba. Depending on the variable, its
1434	!! dimensions are also determined by the number of PFT's (::nvm), number of
1435	!! soil types (::nstm), number of non-vegetative surface types (::nnobio),
1436	!! number of soil levels (::ngrnd), number of soil layers in the hydrological
1437	!! model (i.e. cwrr) (::nslm). Values for these variables are set in
1438	!! constantes_soil.f90 and constantes_veg.f90.\n
1439	!!
1440	!! Memory is allocated for all Sechiba variables and new indexing tables
1441	!! are build making use of both (::kjpij) and (::kjpindex). New indexing tables
1442	!! are needed because a single pixel can contain several PFTs, soil types, etc.
1443	!! The new indexing tables have separate indices for the different
1444	!! PFTs, soil types, etc.\n
1445	!!
1446	!! RECENT CHANGE(S): None
1447	!!
1448	!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
1449	!! variables. However, the routine allocates memory and builds new indexing
1450	!! variables for later use.
1451	!!
1452	!! REFERENCE(S) : None
1453	!!
1454	!! FLOWCHART : None
1455	!! \n
1456	!_ ================================================================================================================================
1457
1458	SUBROUTINE sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)
1459
1460	!! 0.1 Input variables
1461
1462	INTEGER(i_std), INTENT (in) :: kjit !! Time step number (unitless)
1463	INTEGER(i_std), INTENT (in) :: kjpij !! Total size of the un-compressed grid (unitless)
1464	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
1465	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
1466	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless)
1467	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates (latitude,longitude)
1468	!! for pixels (degrees)
1469	!! 0.2 Output variables
1470
1471	!! 0.3 Modified variables
1472
1473	!! 0.4 Local variables
1474
1475	INTEGER(i_std) :: ier !! Check errors in memory allocation (unitless)
1476	INTEGER(i_std) :: ji, jv !! Indeces (unitless)
1477	!_ ==============================================================================================================================
1478
1479	!! 1. Initialize variables
1480
1481	! Dynamic allocation with user-specified dimensions on first call
1482	IF (l_first_sechiba) THEN
1483	l_first_sechiba=.FALSE.
1484	ELSE
1485	CALL ipslerr_p(3,'sechiba_init',' l_first_sechiba false . we stop ','','')
1486	ENDIF
1487
1488	!! Initialize local printlev
1489	printlev_loc=get_printlev('sechiba')
1490
1491
1492	!! 1.1 Initialize 3D vegetation indexation table
1493	ALLOCATE (indexveg(kjpindex*nvm),stat=ier)
1494	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexveg','','')
1495
1496	ALLOCATE (indexlai(kjpindex*(nlai+1)),stat=ier)
1497	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlai','','')
1498
1499	ALLOCATE (indexsoil(kjpindex*nstm),stat=ier)
1500	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsoil','','')
1501
1502	ALLOCATE (indexnobio(kjpindex*nnobio),stat=ier)
1503	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnobio','','')
1504
1505	ALLOCATE (indexgrnd(kjpindex*ngrnd),stat=ier)
1506	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexgrnd','','')
1507
1508	ALLOCATE (indexsnow(kjpindex*nsnow),stat=ier)
1509	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsnow','','')
1510
1511	ALLOCATE (indexlayer(kjpindex*nslm),stat=ier)
1512	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlayer','','')
1513
1514	ALLOCATE (indexnslm(kjpindex*nslm),stat=ier)
1515	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnslm','','')
1516
1517	ALLOCATE (indexalb(kjpindex*2),stat=ier)
1518	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexalb','','')
1519
1520	!! 1.2 Initialize 1D array allocation with restartable value
1521	ALLOCATE (flood_res(kjpindex),stat=ier)
1522	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for flood_res','','')
1523	flood_res(:) = undef_sechiba
1524
1525	ALLOCATE (flood_frac(kjpindex),stat=ier)
1526	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for kjpindex','','')
1527	flood_frac(:) = undef_sechiba
1528
1529	ALLOCATE (snow(kjpindex),stat=ier)
1530	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow','','')
1531	snow(:) = undef_sechiba
1532
1533	ALLOCATE (snow_age(kjpindex),stat=ier)
1534	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_age','','')
1535	snow_age(:) = undef_sechiba
1536
1537	ALLOCATE (drysoil_frac(kjpindex),stat=ier)
1538	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drysoil_frac','','')
1539
1540	ALLOCATE (evap_bare_lim(kjpindex),stat=ier)
1541	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim','','')
1542
1543	ALLOCATE (evap_bare_lim_ns(kjpindex,nstm),stat=ier)
1544	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim_ns','','')
1545
1546	ALLOCATE (evapot(kjpindex),stat=ier)
1547	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot','','')
1548	evapot(:) = undef_sechiba
1549
1550	ALLOCATE (evapot_corr(kjpindex),stat=ier)
1551	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot_corr','','')
1552
1553	ALLOCATE (humrel(kjpindex,nvm),stat=ier)
1554	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for humrel','','')
1555	humrel(:,:) = undef_sechiba
1556
1557	ALLOCATE (vegstress(kjpindex,nvm),stat=ier)
1558	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vegstress','','')
1559	vegstress(:,:) = undef_sechiba
1560
1561	ALLOCATE (njsc(kjpindex),stat=ier)
1562	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for njsc','','')
1563	njsc(:)= undef_int
1564
1565	ALLOCATE (soiltile(kjpindex,nstm),stat=ier)
1566	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soiltile','','')
1567
1568	ALLOCATE (fraclut(kjpindex,nlut),stat=ier)
1569	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fraclut','','')
1570
1571	ALLOCATE (nwdFraclut(kjpindex,nlut),stat=ier)
1572	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for nwdFraclut','','')
1573
1574	ALLOCATE (reinf_slope(kjpindex),stat=ier)
1575	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope','','')
1576
1577	ALLOCATE (reinf_slope_soil(kjpindex, nstm),stat=ier)
1578	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope_soil','','') !
1579
1580	!salma: adding soil params
1581	ALLOCATE (ks(kjpindex),stat=ier)
1582	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ks','','')
1583
1584	ALLOCATE (nvan(kjpindex),stat=ier)
1585	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for nvan ','','')
1586
1587	ALLOCATE (avan(kjpindex),stat=ier)
1588	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for avan','','')
1589
1590	ALLOCATE (mcr(kjpindex),stat=ier)
1591	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcr','','')
1592
1593	ALLOCATE (mcs(kjpindex),stat=ier)
1594	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcs','','')
1595
1596	ALLOCATE (mcfc(kjpindex),stat=ier)
1597	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcfc','','')
1598
1599	ALLOCATE (mcw(kjpindex),stat=ier)
1600	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcw','','')
1601	!end salma: adding soil params
1602
1603
1604
1605
1606	ALLOCATE (vbeta1(kjpindex),stat=ier)
1607	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta1','','')
1608
1609	ALLOCATE (vbeta4(kjpindex),stat=ier)
1610	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta4','','')
1611
1612	ALLOCATE (vbeta5(kjpindex),stat=ier)
1613	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta5','','')
1614
1615	ALLOCATE (soilcap(kjpindex),stat=ier)
1616	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilcap','','')
1617
1618	ALLOCATE (soilflx(kjpindex),stat=ier)
1619	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilflx','','')
1620
1621	ALLOCATE (temp_sol(kjpindex),stat=ier)
1622	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol','','')
1623	temp_sol(:) = undef_sechiba
1624
1625	ALLOCATE (qsurf(kjpindex),stat=ier)
1626	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsurf','','')
1627	qsurf(:) = undef_sechiba
1628
1629	!! 1.3 Initialize 2D array allocation with restartable value
1630	ALLOCATE (qsintveg(kjpindex,nvm),stat=ier)
1631	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintveg','','')
1632	qsintveg(:,:) = undef_sechiba
1633
1634	ALLOCATE (vbeta2(kjpindex,nvm),stat=ier)
1635	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta2','','')
1636
1637	ALLOCATE (vbeta3(kjpindex,nvm),stat=ier)
1638	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3','','')
1639
1640	ALLOCATE (vbeta3pot(kjpindex,nvm),stat=ier)
1641	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3pot','','')
1642
1643	ALLOCATE (gsmean(kjpindex,nvm),stat=ier)
1644	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gsmean','','')
1645
1646	ALLOCATE (cimean(kjpindex,nvm),stat=ier)
1647	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cimean','','')
1648
1649	ALLOCATE (gpp(kjpindex,nvm),stat=ier)
1650	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gpp','','')
1651	gpp(:,:) = undef_sechiba
1652
1653
1654	ALLOCATE (temp_growth(kjpindex),stat=ier)
1655	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_growth','','')
1656	temp_growth(:) = undef_sechiba
1657
1658	ALLOCATE (veget(kjpindex,nvm),stat=ier)
1659	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget','','')
1660	veget(:,:)=undef_sechiba
1661
1662	ALLOCATE (veget_max(kjpindex,nvm),stat=ier)
1663	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget_max','','')
1664
1665	ALLOCATE (tot_bare_soil(kjpindex),stat=ier)
1666	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_bare_soil','','')
1667
1668	ALLOCATE (lai(kjpindex,nvm),stat=ier)
1669	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lai','','')
1670	lai(:,:)=undef_sechiba
1671
1672	ALLOCATE (frac_age(kjpindex,nvm,nleafages),stat=ier)
1673	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_age','','')
1674	frac_age(:,:,:)=undef_sechiba
1675
1676	ALLOCATE (height(kjpindex,nvm),stat=ier)
1677	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for height','','')
1678	height(:,:)=undef_sechiba
1679
1680	ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
1681	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_nobio','','')
1682	frac_nobio(:,:) = undef_sechiba
1683
1684	ALLOCATE (snow_nobio(kjpindex,nnobio),stat=ier)
1685	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio','','')
1686	snow_nobio(:,:) = undef_sechiba
1687
1688	ALLOCATE (snow_nobio_age(kjpindex,nnobio),stat=ier)
1689	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio_age','','')
1690	snow_nobio_age(:,:) = undef_sechiba
1691
1692	ALLOCATE (assim_param(kjpindex,nvm,npco2),stat=ier)
1693	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for assim_param','','')
1694
1695	!! 1.4 Initialize 1D array allocation
1696	ALLOCATE (vevapflo(kjpindex),stat=ier)
1697	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapflo','','')
1698	vevapflo(:)=zero
1699
1700	ALLOCATE (vevapsno(kjpindex),stat=ier)
1701	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapsno','','')
1702
1703	ALLOCATE (vevapnu(kjpindex),stat=ier)
1704	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapnu','','')
1705
1706	ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
1707	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for totfrac_nobio','','')
1708
1709	ALLOCATE (floodout(kjpindex),stat=ier)
1710	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for floodout','','')
1711
1712	ALLOCATE (runoff(kjpindex),stat=ier)
1713	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for runoff','','')
1714
1715	ALLOCATE (drainage(kjpindex),stat=ier)
1716	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drainage','','')
1717
1718	ALLOCATE (returnflow(kjpindex),stat=ier)
1719	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for returnflow','','')
1720	returnflow(:) = zero
1721
1722	ALLOCATE (reinfiltration(kjpindex),stat=ier)
1723	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinfiltration','','')
1724	reinfiltration(:) = zero
1725
1726	ALLOCATE (irrigation(kjpindex),stat=ier)
1727	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for irrigation','','')
1728	irrigation(:) = zero
1729
1730	ALLOCATE (z0h(kjpindex),stat=ier)
1731	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0h','','')
1732
1733	ALLOCATE (z0m(kjpindex),stat=ier)
1734	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0m','','')
1735
1736	ALLOCATE (roughheight(kjpindex),stat=ier)
1737	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for roughheight','','')
1738
1739	ALLOCATE (emis(kjpindex),stat=ier)
1740	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for emis','','')
1741
1742	ALLOCATE (tot_melt(kjpindex),stat=ier)
1743	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_melt','','')
1744
1745	ALLOCATE (vbeta(kjpindex),stat=ier)
1746	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta','','')
1747
1748	ALLOCATE (rau(kjpindex),stat=ier)
1749	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rau','','')
1750
1751	ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
1752	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for deadleaf_cover','','')
1753
1754	ALLOCATE (stempdiag(kjpindex, nslm),stat=ier)
1755	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for stempdiag','','')
1756
1757	ALLOCATE (ftempdiag(kjpindex, ngrnd),stat=ier)
1758	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ftempdiag','','')
1759
1760	ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
1761	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for co2_flux','','')
1762	co2_flux(:,:)=zero
1763
1764	ALLOCATE (shumdiag(kjpindex,nslm),stat=ier)
1765	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag','','')
1766
1767	ALLOCATE (shumdiag_perma(kjpindex,nslm),stat=ier)
1768	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag_perma','','')
1769
1770	ALLOCATE (litterhumdiag(kjpindex),stat=ier)
1771	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for litterhumdiag','','')
1772
1773	ALLOCATE (ptnlev1(kjpindex),stat=ier)
1774	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ptnlev1','','')
1775
1776	ALLOCATE (k_litt(kjpindex),stat=ier)
1777	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for k_litt','','')
1778
1779	!! 1.5 Initialize 2D array allocation
1780	ALLOCATE (vevapwet(kjpindex,nvm),stat=ier)
1781	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapwet','','')
1782	vevapwet(:,:)=undef_sechiba
1783
1784	ALLOCATE (transpir(kjpindex,nvm),stat=ier)
1785	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpir','','')
1786
1787	ALLOCATE (transpot(kjpindex,nvm),stat=ier)
1788	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpot','','')
1789
1790	ALLOCATE (qsintmax(kjpindex,nvm),stat=ier)
1791	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintmax','','')
1792
1793	ALLOCATE (rveget(kjpindex,nvm),stat=ier)
1794	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rveget','','')
1795
1796	ALLOCATE (rstruct(kjpindex,nvm),stat=ier)
1797	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rstruct','','')
1798
1799	ALLOCATE (pgflux(kjpindex),stat=ier)
1800	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for pgflux','','')
1801	pgflux(:)= 0.0
1802
1803	ALLOCATE (cgrnd_snow(kjpindex,nsnow),stat=ier)
1804	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cgrnd_snow','','')
1805	cgrnd_snow(:,:) = 0
1806
1807	ALLOCATE (dgrnd_snow(kjpindex,nsnow),stat=ier)
1808	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for dgrnd_snow','','')
1809	dgrnd_snow(:,:) = 0
1810
1811	ALLOCATE (lambda_snow(kjpindex),stat=ier)
1812	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lambda_snow','','')
1813	lambda_snow(:) = 0
1814
1815	ALLOCATE (temp_sol_add(kjpindex),stat=ier)
1816	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol_add','','')
1817
1818	ALLOCATE (gtemp(kjpindex),stat=ier)
1819	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gtemp','','')
1820
1821	ALLOCATE (frac_snow_veg(kjpindex),stat=ier)
1822	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_snow_veg','','')
1823
1824	ALLOCATE (frac_snow_nobio(kjpindex,nnobio),stat=ier)
1825	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_snow_nobio','','')
1826
1827	ALLOCATE (snowrho(kjpindex,nsnow),stat=ier)
1828	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowrho','','')
1829
1830	ALLOCATE (snowheat(kjpindex,nsnow),stat=ier)
1831	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowheat','','')
1832
1833	ALLOCATE (snowgrain(kjpindex,nsnow),stat=ier)
1834	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowgrain','','')
1835
1836	ALLOCATE (snowtemp(kjpindex,nsnow),stat=ier)
1837	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowtemp','','')
1838
1839	ALLOCATE (snowdz(kjpindex,nsnow),stat=ier)
1840	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowdz','','')
1841
1842	ALLOCATE (mc_layh(kjpindex, nslm),stat=ier)
1843	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mc_layh','','')
1844
1845	ALLOCATE (mcl_layh(kjpindex, nslm),stat=ier)
1846	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcl_layh','','')
1847
1848	ALLOCATE (soilmoist(kjpindex, nslm),stat=ier)
1849	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilmoist','','')
1850
1851
1852	!1.5 Irrigation related variables
1853	ALLOCATE (root_deficit(kjpindex),stat=ier)
1854	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for root_deficit','','') !
1855
1856	ALLOCATE (irrig_frac(kjpindex),stat=ier)
1857	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for irrig_frac','','') !
1858	irrigation(:) = zero
1859
1860	ALLOCATE (irrigated_next(kjpindex),stat=ier) !
1861	IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable irrigated_next','','') !
1862
1863	ALLOCATE (fraction_aeirrig_sw(kjpindex),stat=ier) !
1864	IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fraction_aeirrig_sw','','')
1865
1866	!! 1.6 Initialize indexing table for the vegetation fields.
1867	! In SECHIBA we work on reduced grids but to store in the full 3D filed vegetation variable
1868	! we need another index table : indexveg, indexsoil, indexnobio and indexgrnd
1869	DO ji = 1, kjpindex
1870	!
1871	DO jv = 1, nlai+1
1872	indexlai((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1873	ENDDO
1874	!
1875	DO jv = 1, nvm
1876	indexveg((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1877	ENDDO
1878	!
1879	DO jv = 1, nstm
1880	indexsoil((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1881	ENDDO
1882	!
1883	DO jv = 1, nnobio
1884	indexnobio((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1885	ENDDO
1886	!
1887	DO jv = 1, ngrnd
1888	indexgrnd((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1889	ENDDO
1890	!
1891	DO jv = 1, nsnow
1892	indexsnow((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
1893	ENDDO
1894
1895	DO jv = 1, nslm
1896	indexnslm((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
1897	ENDDO
1898
1899	DO jv = 1, nslm
1900	indexlayer((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1901	ENDDO
1902	!
1903	DO jv = 1, 2
1904	indexalb((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij + offset_omp - offset_mpi
1905	ENDDO
1906	!
1907	ENDDO
1908
1909	!! 2. Read the default value that will be put into variable which are not in the restart file
1910	CALL ioget_expval(val_exp)
1911
1912	IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_init done '
1913
1914	END SUBROUTINE sechiba_init
1915
1916
1917	!! ==============================================================================================================================\n
1918	!! SUBROUTINE : sechiba_clear
1919	!!
1920	!>\BRIEF Deallocate memory of sechiba's variables
1921	!!
1922	!! DESCRIPTION : None
1923	!!
1924	!! RECENT CHANGE(S): None
1925	!!
1926	!! MAIN OUTPUT VARIABLE(S): None
1927	!!
1928	!! REFERENCE(S) : None
1929	!!
1930	!! FLOWCHART : None
1931	!! \n
1932	!_ ================================================================================================================================
1933
1934	SUBROUTINE sechiba_clear()
1935
1936	!! 1. Initialize first run
1937
1938	l_first_sechiba=.TRUE.
1939
1940	!! 2. Deallocate dynamic variables of sechiba
1941
1942	IF ( ALLOCATED (indexveg)) DEALLOCATE (indexveg)
1943	IF ( ALLOCATED (indexlai)) DEALLOCATE (indexlai)
1944	IF ( ALLOCATED (indexsoil)) DEALLOCATE (indexsoil)
1945	IF ( ALLOCATED (indexnobio)) DEALLOCATE (indexnobio)
1946	IF ( ALLOCATED (indexsnow)) DEALLOCATE (indexsnow)
1947	IF ( ALLOCATED (indexgrnd)) DEALLOCATE (indexgrnd)
1948	IF ( ALLOCATED (indexlayer)) DEALLOCATE (indexlayer)
1949	IF ( ALLOCATED (indexnslm)) DEALLOCATE (indexnslm)
1950	IF ( ALLOCATED (indexalb)) DEALLOCATE (indexalb)
1951	IF ( ALLOCATED (flood_res)) DEALLOCATE (flood_res)
1952	IF ( ALLOCATED (flood_frac)) DEALLOCATE (flood_frac)
1953	IF ( ALLOCATED (snow)) DEALLOCATE (snow)
1954	IF ( ALLOCATED (snow_age)) DEALLOCATE (snow_age)
1955	IF ( ALLOCATED (drysoil_frac)) DEALLOCATE (drysoil_frac)
1956	IF ( ALLOCATED (evap_bare_lim)) DEALLOCATE (evap_bare_lim)
1957	IF ( ALLOCATED (evap_bare_lim_ns)) DEALLOCATE (evap_bare_lim_ns)
1958	IF ( ALLOCATED (evapot)) DEALLOCATE (evapot)
1959	IF ( ALLOCATED (evapot_corr)) DEALLOCATE (evapot_corr)
1960	IF ( ALLOCATED (humrel)) DEALLOCATE (humrel)
1961	IF ( ALLOCATED (vegstress)) DEALLOCATE (vegstress)
1962	IF ( ALLOCATED (soiltile)) DEALLOCATE (soiltile)
1963	IF ( ALLOCATED (fraclut)) DEALLOCATE (fraclut)
1964	IF ( ALLOCATED (nwdFraclut)) DEALLOCATE (nwdFraclut)
1965	IF ( ALLOCATED (njsc)) DEALLOCATE (njsc)
1966	IF ( ALLOCATED (reinf_slope)) DEALLOCATE (reinf_slope)
1967	IF ( ALLOCATED (reinf_slope_soil)) DEALLOCATE (reinf_slope_soil)
1968
1969	!salma: adding soil hydraulic params
1970	IF ( ALLOCATED (ks)) DEALLOCATE (ks)
1971	IF ( ALLOCATED (nvan)) DEALLOCATE (nvan)
1972	IF ( ALLOCATED (avan)) DEALLOCATE (avan)
1973	IF ( ALLOCATED (mcr)) DEALLOCATE (mcr)
1974	IF ( ALLOCATED (mcs)) DEALLOCATE (mcs)
1975	IF ( ALLOCATED (mcfc)) DEALLOCATE (mcfc)
1976	IF ( ALLOCATED (mcw)) DEALLOCATE (mcw)
1977	!end salma: adding soil hydraulic params
1978
1979	IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1)
1980	IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4)
1981	IF ( ALLOCATED (vbeta5)) DEALLOCATE (vbeta5)
1982	IF ( ALLOCATED (soilcap)) DEALLOCATE (soilcap)
1983	IF ( ALLOCATED (soilflx)) DEALLOCATE (soilflx)
1984	IF ( ALLOCATED (temp_sol)) DEALLOCATE (temp_sol)
1985	IF ( ALLOCATED (qsurf)) DEALLOCATE (qsurf)
1986	IF ( ALLOCATED (qsintveg)) DEALLOCATE (qsintveg)
1987	IF ( ALLOCATED (vbeta2)) DEALLOCATE (vbeta2)
1988	IF ( ALLOCATED (vbeta3)) DEALLOCATE (vbeta3)
1989	IF ( ALLOCATED (vbeta3pot)) DEALLOCATE (vbeta3pot)
1990	IF ( ALLOCATED (gsmean)) DEALLOCATE (gsmean)
1991	IF ( ALLOCATED (cimean)) DEALLOCATE (cimean)
1992	IF ( ALLOCATED (gpp)) DEALLOCATE (gpp)
1993	IF ( ALLOCATED (temp_growth)) DEALLOCATE (temp_growth)
1994	IF ( ALLOCATED (veget)) DEALLOCATE (veget)
1995	IF ( ALLOCATED (veget_max)) DEALLOCATE (veget_max)
1996	IF ( ALLOCATED (tot_bare_soil)) DEALLOCATE (tot_bare_soil)
1997	IF ( ALLOCATED (lai)) DEALLOCATE (lai)
1998	IF ( ALLOCATED (frac_age)) DEALLOCATE (frac_age)
1999	IF ( ALLOCATED (height)) DEALLOCATE (height)
2000	IF ( ALLOCATED (roughheight)) DEALLOCATE (roughheight)
2001	IF ( ALLOCATED (frac_nobio)) DEALLOCATE (frac_nobio)
2002	IF ( ALLOCATED (snow_nobio)) DEALLOCATE (snow_nobio)
2003	IF ( ALLOCATED (snow_nobio_age)) DEALLOCATE (snow_nobio_age)
2004	IF ( ALLOCATED (assim_param)) DEALLOCATE (assim_param)
2005	IF ( ALLOCATED (vevapflo)) DEALLOCATE (vevapflo)
2006	IF ( ALLOCATED (vevapsno)) DEALLOCATE (vevapsno)
2007	IF ( ALLOCATED (vevapnu)) DEALLOCATE (vevapnu)
2008	IF ( ALLOCATED (totfrac_nobio)) DEALLOCATE (totfrac_nobio)
2009	IF ( ALLOCATED (floodout)) DEALLOCATE (floodout)
2010	IF ( ALLOCATED (runoff)) DEALLOCATE (runoff)
2011	IF ( ALLOCATED (drainage)) DEALLOCATE (drainage)
2012	IF ( ALLOCATED (reinfiltration)) DEALLOCATE (reinfiltration)
2013	IF ( ALLOCATED (irrigation)) DEALLOCATE (irrigation)
2014	IF ( ALLOCATED (tot_melt)) DEALLOCATE (tot_melt)
2015	IF ( ALLOCATED (vbeta)) DEALLOCATE (vbeta)
2016	IF ( ALLOCATED (rau)) DEALLOCATE (rau)
2017	IF ( ALLOCATED (deadleaf_cover)) DEALLOCATE (deadleaf_cover)
2018	IF ( ALLOCATED (stempdiag)) DEALLOCATE (stempdiag)
2019	IF ( ALLOCATED (ftempdiag)) DEALLOCATE (ftempdiag)
2020	IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
2021	IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
2022	IF ( ALLOCATED (shumdiag_perma)) DEALLOCATE (shumdiag_perma)
2023	IF ( ALLOCATED (litterhumdiag)) DEALLOCATE (litterhumdiag)
2024	IF ( ALLOCATED (ptnlev1)) DEALLOCATE (ptnlev1)
2025	IF ( ALLOCATED (k_litt)) DEALLOCATE (k_litt)
2026	IF ( ALLOCATED (vevapwet)) DEALLOCATE (vevapwet)
2027	IF ( ALLOCATED (transpir)) DEALLOCATE (transpir)
2028	IF ( ALLOCATED (transpot)) DEALLOCATE (transpot)
2029	IF ( ALLOCATED (qsintmax)) DEALLOCATE (qsintmax)
2030	IF ( ALLOCATED (rveget)) DEALLOCATE (rveget)
2031	IF ( ALLOCATED (rstruct)) DEALLOCATE (rstruct)
2032	IF ( ALLOCATED (frac_snow_veg)) DEALLOCATE (frac_snow_veg)
2033	IF ( ALLOCATED (frac_snow_nobio)) DEALLOCATE (frac_snow_nobio)
2034	IF ( ALLOCATED (snowrho)) DEALLOCATE (snowrho)
2035	IF ( ALLOCATED (snowgrain)) DEALLOCATE (snowgrain)
2036	IF ( ALLOCATED (snowtemp)) DEALLOCATE (snowtemp)
2037	IF ( ALLOCATED (snowdz)) DEALLOCATE (snowdz)
2038	IF ( ALLOCATED (snowheat)) DEALLOCATE (snowheat)
2039	IF ( ALLOCATED (cgrnd_snow)) DEALLOCATE (cgrnd_snow)
2040	IF ( ALLOCATED (dgrnd_snow)) DEALLOCATE (dgrnd_snow)
2041	IF ( ALLOCATED (lambda_snow)) DEALLOCATE(lambda_snow)
2042	IF ( ALLOCATED (gtemp)) DEALLOCATE (gtemp)
2043	IF ( ALLOCATED (pgflux)) DEALLOCATE (pgflux)
2044	IF ( ALLOCATED (mc_layh)) DEALLOCATE (mc_layh)
2045	IF ( ALLOCATED (mcl_layh)) DEALLOCATE (mcl_layh)
2046	IF ( ALLOCATED (soilmoist)) DEALLOCATE (soilmoist)
2047	IF ( ALLOCATED (root_deficit)) DEALLOCATE (root_deficit)
2048	IF ( ALLOCATED (irrig_frac)) DEALLOCATE (irrig_frac)
2049	IF ( ALLOCATED (irrigated_next)) DEALLOCATE (irrigated_next)
2050
2051	!! 3. Clear all allocated memory
2052
2053	CALL pft_parameters_clear
2054	CALL slowproc_clear
2055	CALL diffuco_clear
2056	CALL enerbil_clear
2057	CALL hydrol_clear
2058	CALL thermosoil_clear
2059	CALL condveg_clear
2060	CALL routing_wrapper_clear
2061
2062	END SUBROUTINE sechiba_clear
2063
2064
2065	!! ==============================================================================================================================\n
2066	!! SUBROUTINE : sechiba_var_init
2067	!!
2068	!>\BRIEF Calculate air density as a function of air temperature and
2069	!! pressure for each terrestrial pixel.
2070	!!
2071	!! RECENT CHANGE(S): None
2072	!!
2073	!! MAIN OUTPUT VARIABLE(S): air density (::rau, kg m^{-3}).
2074	!!
2075	!! REFERENCE(S) : None
2076	!!
2077	!! FLOWCHART : None
2078	!! \n
2079	!_ ================================================================================================================================
2080
2081	SUBROUTINE sechiba_var_init (kjpindex, rau, pb, temp_air)
2082
2083	!! 0.1 Input variables
2084
2085	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
2086	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
2087	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
2088
2089	!! 0.2 Output variables
2090
2091	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rau !! Air density @tex $(kg m^{-3})$ @endtex
2092
2093	!! 0.3 Modified variables
2094
2095	!! 0.4 Local variables
2096
2097	INTEGER(i_std) :: ji !! Indices (unitless)
2098	!_ ================================================================================================================================
2099
2100	!! 1. Calculate intial air density (::rau)
2101
2102	DO ji = 1,kjpindex
2103	rau(ji) = pa_par_hpa * pb(ji) / (cte_molr*temp_air(ji))
2104	END DO
2105
2106	IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_var_init done '
2107
2108	END SUBROUTINE sechiba_var_init
2109
2110
2111	!! ==============================================================================================================================\n
2112	!! SUBROUTINE : sechiba_end
2113	!!
2114	!>\BRIEF Swap old for newly calculated soil temperature.
2115	!!
2116	!! RECENT CHANGE(S): None
2117	!!
2118	!! MAIN OUTPUT VARIABLE(S): soil temperature (::temp_sol; K)
2119	!!
2120	!! REFERENCE(S) : None
2121	!!
2122	!! FLOWCHART : None
2123	!! \n
2124	!! ================================================================================================================================
2125
2126	SUBROUTINE sechiba_end (kjpindex, temp_sol_new, temp_sol)
2127
2128
2129	!! 0.1 Input variables
2130
2131	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
2132	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature (K)
2133
2134	!! 0.2 Output variables
2135	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol !! Soil temperature (K)
2136
2137	!_ ================================================================================================================================
2138
2139	!! 1. Swap temperature
2140
2141	temp_sol(:) = temp_sol_new(:)
2142
2143	IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_end done '
2144
2145	END SUBROUTINE sechiba_end
2146
2147	!! ==============================================================================================================================\n
2148	!! SUBROUTINE : sechiba_interface_orchidee_inca
2149	!!
2150	!>\BRIEF make the interface between surface and atmospheric chemistry
2151	!!
2152	!! DESCRIPTION : This subroutine is called from INCA, the atmospheric chemistry model. It is used to transfer variables from ORCHIDEE to INCA.
2153	!!
2154	!! RECENT CHANGE(S): move from chemistry module to be more generic (feb - 2017)
2155	!!
2156	!! MAIN OUTPUT VARIABLE(S): emission COV to be transport by orchidee to inca in fields_out array
2157	!!
2158	!! REFERENCE(S) : None
2159	!!
2160	!! FLOWCHART : None
2161	!! \n
2162	!! ================================================================================================================================
2163	SUBROUTINE sechiba_interface_orchidee_inca( &
2164	nvm_out, veget_max_out, veget_frac_out, lai_out, snow_out, &
2165	field_out_names, fields_out, field_in_names, fields_in)
2166
2167
2168	INTEGER, INTENT(out) :: nvm_out !! Number of vegetation types
2169	REAL(r_std), DIMENSION (:,:), INTENT(out) :: veget_max_out !! Max. fraction of vegetation type (LAI -> infty)
2170	REAL(r_std), DIMENSION (:,:), INTENT(out) :: veget_frac_out !! Fraction of vegetation type (unitless, 0-1)
2171	REAL(r_std), DIMENSION (:,:), INTENT(out) :: lai_out !! Surface foliere
2172	REAL(r_std), DIMENSION (:) , INTENT(out) :: snow_out !! Snow mass [Kg/m^2]
2173
2174	!
2175	! Optional arguments
2176	!
2177	! Names and fields for emission variables : to be transport by Orchidee to Inca
2178	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_out_names
2179	REAL(r_std),DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: fields_out
2180	!
2181	! Names and fields for deposit variables : to be transport from chemistry model by INCA to ORCHIDEE.
2182	CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_in_names
2183	REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(IN) :: fields_in
2184
2185
2186	! Variables always transmitted from sechiba to inca
2187	nvm_out = nvm
2188	veget_max_out(:,:) = veget_max(:,:)
2189	veget_frac_out(:,:) = veget(:,:)
2190	lai_out(:,:) = lai(:,:)
2191	snow_out(:) = snow(:)
2192
2193	! Call chemistry_flux_interface if at least one of variables field_out_names or
2194	! field_in_names is present in the argument list of sechiba_interface_orchidee_inca when called from inca.
2195	IF (PRESENT(field_out_names) .AND. .NOT. PRESENT(field_in_names)) THEN
2196	CALL chemistry_flux_interface(field_out_names=field_out_names, fields_out=fields_out)
2197	ELSE IF (.NOT. PRESENT(field_out_names) .AND. PRESENT(field_in_names)) THEN
2198	CALL chemistry_flux_interface(field_in_names=field_in_names, fields_in=fields_in)
2199	ELSE IF (PRESENT(field_out_names) .AND. PRESENT(field_in_names)) THEN
2200	CALL chemistry_flux_interface(field_out_names=field_out_names, fields_out=fields_out, &
2201	field_in_names=field_in_names, fields_in=fields_in)
2202	ENDIF
2203
2204	END SUBROUTINE sechiba_interface_orchidee_inca
2205
2206
2207	END MODULE sechiba
2208

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

Download in other formats: