Context Navigation

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

Last change on this file since 7709 was 7709, checked in by josefine.ghattas, 2 years ago

Integrated new irrigation scheme developed by Pedro Arboleda. See ticket #857
This corresponds to revsion 7708 of version pedro.arboleda/ORCHIDEE. Following differences were made but were not made on the pedro.arboleda/ORCHIDEE :

argumet place in call to routing_wrapper_intialize changed order
lines with only change in space were not taken
some indentation changed
set irrigation output as enalbled false if not do_irrigation

Property svn:keywords set to Revision Date HeadURL Date Author Revision

File size: 116.5 KB

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

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

Download in other formats: