Context Navigation

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

Last change on this file since 5820 was 5820, checked in by josefine.ghattas, 5 years ago
sechiba.f90: Cut line to long for gfortran makeorchidee_fcm: small change to remove error message when compiling first time
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 111.6 KB

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

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

Download in other formats: