Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sechiba.f90 @ 5501

Last change on this file since 5501 was 5501, checked in by aclsce, 4 years ago
First import of IPSLCM6.5_work_ENSEMBLES working configuration
File size: 111.5 KB

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

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

Context Navigation

source: CONFIG_DEVT/IPSLCM6.5_work_ENSEMBLES/modeles/ORCHIDEE/src_sechiba/sechiba.f90 @ 5501

Download in other formats: