Context Navigation

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

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

Correction to avoid problems with negative soil moisture. See ticket #325

Ducharne

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

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

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

Download in other formats: