Context Navigation

source: tags/ORCHIDEE_2_1/ORCHIDEE/src_sechiba/sechiba.f90 @ 5630

Last change on this file since 5630 was 5573, checked in by josefine.ghattas, 6 years ago
Do not use t2m/q2m coming from the atmospheric model anymore and instead use temp_air and qair from first atmpospheric model layer. t2m/q2m were previously used only in diffuco_trans_co2 and in stomate_initialize. Also removed diagnostic variables t2m and q2m. See ticket #372
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 111.2 KB

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

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

Download in other formats: