Context Navigation

sechiba.f90 @ 7746

Last change on this file since 7746 was 4074, checked in by jan.polcher, 8 years ago
Convergence with Trunk version 4061 and in particular introduces the option FROZ_FRAC_CORR to reduce runoff over frozen soils.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 106.7 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_sechiba/sechiba.f90 @ 7746

Download in other formats: