source: branches/publications/ORCHIDEE-MICT-OP-r6850/src_sechiba/sechiba.f90 @ 8005

Last change on this file since 8005 was 6849, checked in by yidi.xu, 4 years ago

ORCHIDEE-MICT-OP for oil palm growth modelling

  • Property svn:keywords set to Revision Date HeadURL Date Author Revision
File size: 161.9 KB
Line 
1
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, hydrolc_main (or hydrol_main),
11!! enerbil_fusion, condveg_main and thermosoil_main. Note that sechiba_main
12!! calls slowproc_main and thus indirectly calculates the biogeochemical
13!! processes as well.
14!!
15!!\n DESCRIPTION  : :: shumdiag, :: litterhumdiag and :: stempdiag are not
16!! saved in the restart file because at the first time step because they
17!! are recalculated. However, they must be saved as they are in slowproc
18!! which is called before the modules which calculate them.
19!!
20!! RECENT CHANGE(S): None
21!!
22!! REFERENCE(S) : None
23!!   
24!! SVN     :
25!! $HeadURL$
26!! $Date$
27!! $Revision$
28!! \n
29!_ ================================================================================================================================
30 
31MODULE 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 hydrolc
45  USE thermosoil
46  USE thermosoilc
47  USE sechiba_io_p
48  USE slowproc
49  USE routing
50  USE ioipsl_para
51  USE chemistry
52  USE vegetation
53
54  IMPLICIT NONE
55
56  PRIVATE
57  PUBLIC sechiba_main, sechiba_initialize, sechiba_clear, sechiba_interface_orchidee_inca
58
59  INTEGER(i_std), SAVE                             :: printlev_loc   !! local printlev for this module
60!$OMP THREADPRIVATE(printlev_loc)
61  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexveg       !! indexing array for the 3D fields of vegetation
62!$OMP THREADPRIVATE(indexveg)
63  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlai       !! indexing array for the 3D fields of vegetation for nlai+1
64!$OMP THREADPRIVATE(indexlai)
65  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlai0      !! indexing array for the 3D fields of vegetation for nlai
66!$OMP THREADPRIVATE(indexlai0)
67  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnobio     !! indexing array for the 3D fields of other surfaces (ice,
68                                                                     !! lakes, ...)
69!$OMP THREADPRIVATE(indexnobio)
70  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsoil      !! indexing array for the 3D fields of soil types (kjpindex*nstm)
71!$OMP THREADPRIVATE(indexsoil)
72  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexgrnd      !! indexing array for the 3D ground heat profiles (kjpindex*ngrnd)
73!$OMP THREADPRIVATE(indexgrnd)
74  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlayer     !! indexing array for the 3D fields of soil layers in CWRR (kjpindex*nslm)
75!$OMP THREADPRIVATE(indexlayer)
76  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnslm      !! indexing array for the 3D fields of diagnostic soil layers (kjpindex*nslm)
77!$OMP THREADPRIVATE(indexnslm)
78  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexalb       !! indexing array for the 2 fields of albedo
79!$OMP THREADPRIVATE(indexalb)
80  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsnow      !! indexing array for the 3D fields snow layers
81!$OMP THREADPRIVATE(indexsnow)
82  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: veget          !! Fraction of vegetation type (unitless, 0-1)       
83!$OMP THREADPRIVATE(veget)
84  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: veget_max      !! Max. fraction of vegetation type (LAI -> infty, unitless)
85!$OMP THREADPRIVATE(veget_max)
86  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: tot_bare_soil  !! Total evaporating bare soil fraction
87!$OMP THREADPRIVATE(tot_bare_soil)
88  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: height         !! Vegetation Height (m)
89!$OMP THREADPRIVATE(height)
90  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: totfrac_nobio  !! Total fraction of continental ice+lakes+cities+...
91                                                                     !! (unitless, 0-1)
92!$OMP THREADPRIVATE(totfrac_nobio)
93  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: floodout       !! Flow out of floodplains from hydrol
94!$OMP THREADPRIVATE(floodout)
95  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: runoff         !! Surface runoff calculated by hydrol or hydrolc
96                                                                     !! @tex $(kg m^{-2})$ @endtex
97!$OMP THREADPRIVATE(runoff)
98  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: drainage       !! Deep drainage calculatedd by hydrol or hydrolc
99                                                                     !! @tex $(kg m^{-2})$ @endtex
100!$OMP THREADPRIVATE(drainage)
101  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: returnflow     !! Water flow from lakes and swamps which returns to
102                                                                     !! the grid box @tex $(kg m^{-2})$ @endtex
103!$OMP THREADPRIVATE(returnflow)
104  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: reinfiltration !! Routed water which returns into the soil
105!$OMP THREADPRIVATE(reinfiltration)
106  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: irrigation     !! Irrigation flux taken from the routing reservoirs and
107                                                                     !! being put into the upper layers of the soil
108                                                                     !! @tex $(kg m^{-2})$ @endtex
109!$OMP THREADPRIVATE(irrigation)
110  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: emis           !! Surface emissivity (unitless)
111!$OMP THREADPRIVATE(emis)
112  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: z0h           !! Surface roughness for heat (m)
113!$OMP THREADPRIVATE(z0h)
114  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: z0m           !! Surface roughness for momentum (m)
115!$OMP THREADPRIVATE(z0m)
116  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: roughheight    !! Effective height for roughness (m)
117!$OMP THREADPRIVATE(roughheight)
118  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: roughheight_pft    !! Effective height for roughness (m)
119!$OMP THREADPRIVATE(roughheight_pft)
120  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: reinf_slope    !! slope coefficient (reinfiltration)
121!$OMP THREADPRIVATE(reinf_slope)
122  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: shumdiag       !! Mean relative soil moisture in the different levels used
123                                                                     !! by thermosoil.f90 (unitless, 0-1)
124!$OMP THREADPRIVATE(shumdiag)
125  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: shumdiag_perma !! Saturation degree of the soil
126!$OMP THREADPRIVATE(shumdiag_perma)
127  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: k_litt         !! litter cond.
128!$OMP THREADPRIVATE(k_litt)
129  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: litterhumdiag  !! Litter dryness factor (unitless, 0-1)
130!$OMP THREADPRIVATE(litterhumdiag)
131  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: stempdiag      !! Temperature which controls canopy evolution (K)
132!$OMP THREADPRIVATE(stempdiag)
133  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: qsintveg       !! Water on vegetation due to interception
134                                                                     !! @tex $(kg m^{-2})$ @endtex
135!$OMP THREADPRIVATE(qsintveg)
136  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: tq_cdrag_pft   !! Interception resistance (unitless,0-1)
137!$OMP THREADPRIVATE(tq_cdrag_pft)
138  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta2         !! Interception resistance (unitless,0-1)
139!$OMP THREADPRIVATE(vbeta2)
140  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta3         !! Vegetation resistance (unitless,0-1)
141!$OMP THREADPRIVATE(vbeta3)
142  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta3pot      !! Potential vegetation resistance
143!$OMP THREADPRIVATE(vbeta3pot)
144  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: gsmean         !! Mean stomatal conductance for CO2 (mol m-2 s-1)
145!$OMP THREADPRIVATE(gsmean)
146  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: cimean         !! STOMATE: mean intercellular CO2 concentration (ppm)
147!$OMP THREADPRIVATE(cimean)
148  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vevapwet       !! Interception loss over each PFT
149                                                                     !! @tex $(kg m^{-2} days^{-1})$ @endtex
150!$OMP THREADPRIVATE(vevapwet)
151  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: transpir       !! Transpiration @tex $(kg m^{-2} days^{-1})$ @endtex
152!$OMP THREADPRIVATE(transpir)
153  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: transpot       !! Potential Transpiration (needed for irrigation)
154!$OMP THREADPRIVATE(transpot)
155  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: qsintmax       !! Maximum amount of water in the canopy interception
156                                                                     !! reservoir @tex $(kg m^{-2})$ @endtex
157!$OMP THREADPRIVATE(qsintmax)
158  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: rveget         !! Surface resistance for the vegetation
159                                                                     !! @tex $(s m^{-1})$ @endtex
160!$OMP THREADPRIVATE(rveget)
161  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: rstruct        !! Vegetation structural resistance
162!$OMP THREADPRIVATE(rstruct)
163  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: snow_nobio     !! Snow mass of non-vegetative surfaces
164                                                                     !! @tex $(kg m^{-2})$ @endtex
165!$OMP THREADPRIVATE(snow_nobio)
166  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: snow_nobio_age !! Snow age on non-vegetative surfaces (days)
167!$OMP THREADPRIVATE(snow_nobio_age)
168  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: frac_nobio     !! Fraction of non-vegetative surfaces (continental ice,
169                                                                     !! lakes, ...) (unitless, 0-1)
170!$OMP THREADPRIVATE(frac_nobio)
171  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: assim_param    !! min+max+opt temps, vcmax, vjmax for photosynthesis
172!$OMP THREADPRIVATE(assim_param)
173  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lai            !! Surface foliaire
174!@chy it should be declared here.@xchen
175!$OMP THREADPRIVATE(lai)
176  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lai_cl1            !! @xchen Surface foliaire
177!$OMP THREADPRIVATE(lai_cl1)
178  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lai_cl2            !! @xchen Surface foliaire
179!$OMP THREADPRIVATE(lai_cl2)
180  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lai_cl3            !! @xchen Surface foliaire
181!$OMP THREADPRIVATE(lai_cl3)
182  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lai_cl4            !! @xchen Surface foliaire
183!$OMP THREADPRIVATE(lai_cl4)
184  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: laitop            !! @xchen
185!$OMP THREADPRIVATE(laish)
186  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: laidow            !! @xchen
187!! xchen add vpd for gpp simulation in Amazon
188  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)  :: VPD              !! Vapor Pressure Deficit (kPa)
189  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)  :: ave_VPD              !! Vapor Pressure Deficit (kPa)
190  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vcmax_cl1            !! @xchen Surface foliaire
191  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vcmax_cl2            !! @xchen Surface foliaire
192  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vcmax_cl3            !! @xchen Surface foliaire
193  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vcmax_cl4            !! @xchen Surface foliaire
194
195  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: gpp            !! STOMATE: GPP. gC/m**2 of total area
196!$OMP THREADPRIVATE(gpp)
197  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)     :: temp_growth       !! Growth temperature (°C) - Is equal to t2m_month
198!$OMP THREADPRIVATE(temp_growth)
199  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: humrel         !! Relative humidity
200!$OMP THREADPRIVATE(humrel)
201  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vegstress      !! Vegetation moisture stress (only for vegetation growth)
202!$OMP THREADPRIVATE(vegstress)
203  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vegstress_old  !! Vegetation moisture stress of previous time step
204!$OMP THREADPRIVATE(vegstress_old)
205  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: frac_age       !! Age efficacity from STOMATE for isoprene
206!$OMP THREADPRIVATE(frac_age)
207  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: soiltile       !! Fraction of each soil tile (0-1, unitless)
208!$OMP THREADPRIVATE(soiltile)
209  LOGICAL, ALLOCATABLE, SAVE, DIMENSION (:)        :: is_crop_soil   !! whether the soil tile is under cropland
210!$OMP THREADPRIVATE(is_crop_soil)
211  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: fraclut        !! Fraction of each landuse tile (0-1, unitless)
212!$OMP THREADPRIVATE(fraclut)
213  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: nwdFraclut     !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless)
214!$OMP THREADPRIVATE(nwdFraclut)
215  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: njsc           !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
216!$OMP THREADPRIVATE(njsc)
217  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta1         !! Snow resistance
218!$OMP THREADPRIVATE(vbeta1)
219  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta4         !! Bare soil resistance
220!$OMP THREADPRIVATE(vbeta4)
221  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta4_pft     !! Bare soil resistance
222!$OMP THREADPRIVATE(vbeta4_pft)
223  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta5         !! Floodplains resistance
224!$OMP THREADPRIVATE(vbeta5)
225  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: soilcap        !!
226!$OMP THREADPRIVATE(soilcap)
227  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: soilcap_pft    !!
228!$OMP THREADPRIVATE(soilcap_pft)
229  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: soilflx        !!
230!$OMP THREADPRIVATE(soilflx)
231  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: soilflx_pft    !!
232!$OMP THREADPRIVATE(soilflx_pft)
233  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: temp_sol       !! Soil temperature
234!$OMP THREADPRIVATE(temp_sol)
235  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: temp_sol_pft    !! Soil temperature for each pft
236!$OMP THREADPRIVATE(temp_sol_pft)
237  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: temp_sol_new_pft    !! Soil temperature
238!$OMP THREADPRIVATE(temp_sol_new_pft)
239  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: qsurf          !! near soil air moisture
240!$OMP THREADPRIVATE(qsurf)
241  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: flood_res      !! flood reservoir estimate
242!$OMP THREADPRIVATE(flood_res)
243  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: flood_frac     !! flooded fraction
244!$OMP THREADPRIVATE(flood_frac)
245  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: snow           !! Snow mass [Kg/m^2]
246!$OMP THREADPRIVATE(snow)
247  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: snow_age       !! Snow age
248!$OMP THREADPRIVATE(snow_age)
249  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: drysoil_frac   !! Fraction of visibly (albedo) Dry soil (Between 0 and 1)
250!$OMP THREADPRIVATE(drysoil_frac)
251  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: rsol           !! resistance to bare soil evaporation
252!$OMP THREADPRIVATE(rsol)
253  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evap_bare_lim  !! Bare soil stress
254!$OMP THREADPRIVATE(evap_bare_lim)
255!!  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: evap_bare_lim_pft  !! Bare soil stress
256!!!$OMP THREADPRIVATE(evap_bare_lim_pft)
257  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: soil_deficit   !! water deficit to reach IRRIG_FULFILL
258!$OMP THREADPRIVATE(soil_deficit)
259
260  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: co2_flux       !! CO2 flux (gC/m**2 of average ground/s)
261!$OMP THREADPRIVATE(co2_flux)
262  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evapot         !! Soil Potential Evaporation
263!$OMP THREADPRIVATE(evapot)
264  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evapot_corr    !! Soil Potential Evaporation Correction (Milly 1992)
265!$OMP THREADPRIVATE(evapot_corr)
266  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapflo       !! Floodplains evaporation
267!$OMP THREADPRIVATE(vevapflo)
268  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapsno       !! Snow evaporation
269!$OMP THREADPRIVATE(vevapsno)
270  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapnu        !! Bare soil evaporation
271!$OMP THREADPRIVATE(vevapnu)
272  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vevapnu_pft        !! Bare soil evaporation
273!$OMP THREADPRIVATE(vevapnu_pft)
274  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: tot_melt       !! Total melt
275!$OMP THREADPRIVATE(tot_melt)
276  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta          !! Resistance coefficient
277!$OMP THREADPRIVATE(vbeta)
278  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta_pft      !! Resistance coefficient for each pft
279!$OMP THREADPRIVATE(vbeta_pft)
280  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: fusion         !!
281!$OMP THREADPRIVATE(fusion)
282  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: rau            !! Density
283!$OMP THREADPRIVATE(rau)
284  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: deadleaf_cover !! Fraction of soil covered by dead leaves
285!$OMP THREADPRIVATE(deadleaf_cover)
286  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: ptnlev1        !! 1st level Different levels soil temperature
287!$OMP THREADPRIVATE(ptnlev1)
288  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: irrig_frac        !!irrigated fraction of the croplands
289!$OMP THREADPRIVATE(irrig_frac)
290  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: tot_vegfrac_crop !! Total fraction occcupied by crops (0-1, unitess)
291!$OMP THREADPRIVATE(tot_vegfrac_crop)
292  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: mc_layh        !! Volumetric soil moisture for each layer in hydrol(liquid + ice) (m3/m3)
293!$OMP THREADPRIVATE(mc_layh)
294  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: mc_layh_s    !! Volumetric soil moisture for each layer in hydrol(liquid + ice) (m3/m3)
295!$OMP THREADPRIVATE(mc_layh_s)
296!  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: mc_layh_pft    !! Volumetric soil moisture for each layer in hydrol(liquid + ice) (m3/m3)
297!!$OMP THREADPRIVATE(mc_layh_pft)
298  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: mcl_layh       !! Volumetric soil moisture for each layer in hydrol(liquid) (m3/m3)
299!$OMP THREADPRIVATE(mcl_layh)
300  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mcl_layh_s    !! Volumetric soil moisture for each layer in hydrol(liquid) (m3/m3)
301!$OMP THREADPRIVATE(mcl_layh_s)
302  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: soilmoist      !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
303!$OMP THREADPRIVATE(soilmoist)
304
305  LOGICAL, SAVE                                    :: l_first_sechiba = .TRUE. !! Flag controlling the intialisation (true/false)
306!$OMP THREADPRIVATE(l_first_sechiba)
307
308  ! Variables related to snow processes calculations 
309
310  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: frac_snow_veg   !! Snow cover fraction on vegetation (unitless)
311!$OMP THREADPRIVATE(frac_snow_veg)
312  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: frac_snow_nobio !! Snow cover fraction on continental ice, lakes, etc (unitless)
313!$OMP THREADPRIVATE(frac_snow_nobio)
314  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowrho      !! snow density for each layer
315!$OMP THREADPRIVATE(snowrho)
316  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowheat     !! snow heat content for each layer (J/m2)
317!$OMP THREADPRIVATE(snowheat)
318  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowgrain    !! snow grain size (m)
319!$OMP THREADPRIVATE(snowgrain)
320  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowtemp     !! snow temperature profile (K)
321!$OMP THREADPRIVATE(snowtemp)
322  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowdz       !! snow layer thickness (m)
323!$OMP THREADPRIVATE(snowdz)
324  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: gtemp        !! soil surface temperature
325!$OMP THREADPRIVATE(gtemp)
326  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: soilflxresid !! Net flux to the snowpack
327!$OMP THREADPRIVATE(soilflxresid)
328  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: pgflux       !! net energy into snow pack
329!$OMP THREADPRIVATE(pgflux)
330  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: cgrnd_snow   !! Integration coefficient for snow numerical scheme
331!$OMP THREADPRIVATE(cgrnd_snow)
332  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: dgrnd_snow   !! Integration coefficient for snow numerical scheme
333!$OMP THREADPRIVATE(dgrnd_snow)
334  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: lambda_snow  !! Coefficient of the linear extrapolation of surface temperature
335                                                                  !! from the first and second snow layers
336!$OMP THREADPRIVATE(lambda_snow)
337  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
338!$OMP THREADPRIVATE(temp_sol_add)
339  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: sfcfrz       !! snow surface layer frozen fraction
340!$OMP THREADPRIVATE(sfcfrz)
341  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: radsink      !! radiation sink/source (J/m2)
342!$OMP THREADPRIVATE(radsink)
343
344  ! Variables related to deep permafrost calculations
345  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: tdeep          !! deep temperature profile
346!$OMP THREADPRIVATE(tdeep)
347  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: hsdeep         !! deep soil humidity profile
348!$OMP THREADPRIVATE(hsdeep)
349  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: heat_Zimov     !! heating associated with decomposition
350!$OMP THREADPRIVATE(heat_Zimov)
351  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)     :: sfluxCH4_deep    !! surface flux of CH4 to atmosphere from permafrost
352!$OMP THREADPRIVATE(sfluxCH4_deep)
353  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)     :: sfluxCO2_deep    !! surface flux of CO2 to atmosphere from permafrost
354!$OMP THREADPRIVATE(sfluxCO2_deep)
355  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)     :: thawed_humidity 
356!$OMP THREADPRIVATE(thawed_humidity)
357  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)     :: depth_organic_soil
358!$OMP THREADPRIVATE(depth_organic_soil)
359  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)     :: zz_deep
360!$OMP THREADPRIVATE(zz_deep)
361  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)     :: zz_coef_deep
362!$OMP THREADPRIVATE(zz_coef_deep)
363  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: soilc_total    !! total  soil carbon for use in thermal calcs
364!$OMP THREADPRIVATE(soilc_total)
365  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: rot_cmd
366!$OMP THREADPRIVATE(rot_cmd)
367  LOGICAL, ALLOCATABLE, SAVE, DIMENSION(:)        :: f_rot_sech
368!$OMP THREADPRIVATE(f_rot_sech)
369  ! the rotation command matrix: xxxyyzz, xxx: % land fraction change, yy:
370  ! source PFT, zz: destination PFT,  dimensions: kjpindex, rot_cmd_max, cyc_rot_max
371!gmjc top 5 layer grassland soil moisture for grazing
372  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: tmc_topgrass
373!$OMP THREADPRIVATE(tmc_topgrass)
374!end gmjc
375
376  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: humcste_use
377!$OMP THREADPRIVATE(humcste_use)
378  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: altmax
379!$OMP THREADPRIVATE(altmax)
380
381CONTAINS
382
383!!  =============================================================================================================================
384!! SUBROUTINE:    sechiba_initialize
385!!
386!>\BRIEF          Initialize all prinicipal modules by calling their "_initialize" subroutines
387!!
388!! DESCRIPTION:   Initialize all prinicipal modules by calling their "_initialize" subroutines
389!!
390!! \n
391!_ ==============================================================================================================================
392
393  SUBROUTINE sechiba_initialize( &
394       kjit,         kjpij,        kjpindex,     index,      date0,       &
395       lalo,         contfrac,     neighbours,   resolution, zlev,        &
396       u,            v,            qair,         t2m,        temp_air,    &
397       petAcoef,     peqAcoef,     petBcoef,     peqBcoef,                &
398       precip_rain,  precip_snow,  lwdown,       swnet,      swdown,      &
399       pb,           rest_id,      hist_id,      hist2_id,                &
400       rest_id_stom, hist_id_stom, hist_id_stom_IPCC,                     &
401       coastalflow,  riverflow,    tsol_rad,     vevapp,       qsurf_out, &
402       z0m_out,      z0h_out,      albedo,       fluxsens,     fluxlat,      emis_out,  &
403       netco2flux,   fco2_lu,      temp_sol_new,  tq_cdrag)
404
405!! 0.1 Input variables
406    INTEGER(i_std), INTENT(in)                               :: kjit              !! Time step number (unitless)
407    INTEGER(i_std), INTENT(in)                               :: kjpij             !! Total size of the un-compressed grid
408                                                                                  !! (unitless)
409    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size - terrestrial pixels only
410                                                                                  !! (unitless)
411    INTEGER(i_std),INTENT (in)                               :: rest_id           !! _Restart_ file identifier (unitless)
412    INTEGER(i_std),INTENT (in)                               :: hist_id           !! _History_ file identifier (unitless)
413    INTEGER(i_std),INTENT (in)                               :: hist2_id          !! _History_ file 2 identifier (unitless)
414    INTEGER(i_std),INTENT (in)                               :: rest_id_stom      !! STOMATE's _Restart_ file identifier
415                                                                                  !! (unitless)
416    INTEGER(i_std),INTENT (in)                               :: hist_id_stom      !! STOMATE's _History_ file identifier
417                                                                                  !! (unitless)
418    INTEGER(i_std),INTENT(in)                                :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file
419                                                                                  !! identifier (unitless)
420    REAL(r_std), INTENT (in)                                 :: date0             !! Initial date (??unit??)
421    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)          :: lalo              !! Geographic coordinates (latitude,longitude)
422                                                                                  !! for grid cells (degrees)
423    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: contfrac          !! Fraction of continent in the grid
424                                                                                  !! (unitless, 0-1)
425    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index             !! Indices of the pixels on the map.
426                                                                                  !! Sechiba uses a reduced grid excluding oceans
427                                                                                  !! ::index contains the indices of the
428                                                                                  !! terrestrial pixels only! (unitless)
429    INTEGER(i_std), DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours        !! Neighboring grid points if land!(unitless)
430    REAL(r_std), DIMENSION (kjpindex,2), INTENT(in)          :: resolution        !! Size in x and y of the grid (m)
431   
432    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u                 !! Lowest level wind speed in direction u
433                                                                                  !! @tex $(m.s^{-1})$ @endtex
434    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: v                 !! Lowest level wind speed in direction v
435                                                                                  !! @tex $(m.s^{-1})$ @endtex
436    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: zlev              !! Height of first layer (m)
437    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair              !! Lowest level specific humidity
438                                                                                  !! @tex $(kg kg^{-1})$ @endtex
439    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: t2m               !! 2m air temperature (K)
440    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: precip_rain       !! Rain precipitation
441                                                                                  !! @tex $(kg m^{-2})$ @endtex
442    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: precip_snow       !! Snow precipitation
443                                                                                  !! @tex $(kg m^{-2})$ @endtex
444    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: lwdown            !! Down-welling long-wave flux
445                                                                                  !! @tex $(W m^{-2})$ @endtex
446    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swnet             !! Net surface short-wave flux
447                                                                                  !! @tex $(W m^{-2})$ @endtex
448    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swdown            !! Down-welling surface short-wave flux
449                                                                                  !! @tex $(W m^{-2})$ @endtex
450    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_air          !! Air temperature (K)
451    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petAcoef          !! Coefficients A for T from the Planetary
452                                                                                  !! Boundary Layer
453    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqAcoef          !! Coefficients A for q from the Planetary
454                                                                                  !! Boundary Layer
455    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petBcoef          !! Coefficients B for T from the Planetary
456                                                                                  !! Boundary Layer
457    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqBcoef          !! Coefficients B for q from the Planetary
458                                                                                  !! Boundary Layer
459    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: pb                !! Surface pressure (hPa)
460
461
462!! 0.2 Output variables
463    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: coastalflow       !! Outflow on coastal points by small basins.
464                                                                                  !! This is the water which flows in a disperse
465                                                                                  !! way into the ocean
466                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
467    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: riverflow         !! Outflow of the major rivers.
468                                                                                  !! The flux will be located on the continental
469                                                                                  !! grid but this should be a coastal point 
470                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
471    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: tsol_rad          !! Radiative surface temperature
472                                                                                  !! @tex $(W m^{-2})$ @endtex
473    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapp            !! Total of evaporation
474                                                                                  !! @tex $(kg m^{-2} days^{-1})$ @endtex
475   
476    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qsurf_out         !! Surface specific humidity
477                                                                                  !! @tex $(kg kg^{-1})$ @endtex
478    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0m_out           !! Surface roughness momentum (output diagnostic, m)
479    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0h_out           !! Surface roughness heat (output diagnostic, m)
480    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out)         :: albedo            !! Surface albedo for visible and near-infrared
481                                                                                  !! (unitless, 0-1)
482    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxsens          !! Sensible heat flux
483                                                                                  !! @tex $(W m^{-2})$ @endtex
484    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxlat           !! Latent heat flux
485                                                                                  !! @tex $(W m^{-2})$ @endtex
486    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: emis_out          !! Emissivity (unitless)
487    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: netco2flux        !! Sum CO2 flux over PFTs
488                                                                                  !! ??(gC m^{-2} s^{-1})??
489    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fco2_lu           !! Land Cover Change CO2 flux
490                                                                                  !! ??(gC m^{-2} s^{-1})??
491    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: temp_sol_new      !! New ground temperature (K)
492!    REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (out)      :: temp_sol_new_pft  !! New ground temperature (K)
493
494!! 0.3 Modified
495    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: tq_cdrag          !! Surface drag coefficient (-)
496
497!! 0.4 Local variables
498    INTEGER(i_std)                                           :: ji, jv            !! Index (unitless)
499    INTEGER(i_std)                                           :: ier
500    REAL(r_std), DIMENSION(kjpindex)                         :: histvar           !! Computations for history files (unitless)
501    REAL(r_std), DIMENSION(kjpindex)                         :: zmaxh_glo         !! 2D field of constant soil depth (zmaxh) (m)
502    CHARACTER(LEN=80)                                        :: var_name          !! To store variables names for I/O (unitless)
503
504    REAL(r_std), DIMENSION(kjpindex,nslm,nvm)                :: mc_layh_pft
505    REAL(r_std), DIMENSION(kjpindex,nslm,nvm)                :: mcl_layh_pft
506    REAL(r_std), DIMENSION(kjpindex,nslm,nvm)                :: soilmoist_pft 
507    REAL(r_std), DIMENSION(kjpindex,nvm)                     :: veget_max_bg !! No nobio. All placed in bareground pft
508!_ ================================================================================================================================
509
510    IF (printlev>=3) WRITE(numout,*) 'Start sechiba_initialize'
511
512    !! 1. Initialize variables on first call
513
514    !! 1.1 Initialize most of sechiba's variables
515    CALL sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)
516   
517    !! 1.3 Initialize stomate's variables
518    CALL slowproc_initialize (kjit,     kjpij,          kjpindex,       date0,          &     
519                        index,          indexveg,       lalo,           neighbours,     &
520                        resolution,     contfrac,       soiltile,       reinf_slope,    &
521                        t2m,                                                            &           
522                        deadleaf_cover, assim_param,    lai,                            &
523!!@xchen
524                        lai_cl1,        lai_cl2,        lai_cl3,        lai_cl4, laitop, laidow,   VPD, ave_VPD, &
525                        vcmax_cl1,      vcmax_cl2,      vcmax_cl3,      vcmax_cl4,      & 
526                        frac_age,       &
527                        height,         veget,          frac_nobio,     njsc,           & 
528                        veget_max,      totfrac_nobio,  qsintmax,       rest_id,        &
529                        rest_id_stom,   hist_id_stom,   tot_bare_soil,  fraclut,        &
530                        nwdfraclut,                                                     &
531                        hist_id_stom_IPCC, co2_flux,    fco2_lu,                        &
532                        temp_growth,                                                    &
533                        soilc_total,    thawed_humidity, depth_organic_soil,   heat_Zimov, f_rot_sech,altmax)
534
535    netco2flux(:) = zero
536    DO jv = 2,nvm
537       netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*veget_max(:,jv)
538    ENDDO
539   
540    !!! xuhui: allocation of rot_cmd has to be after slowproc_initialize
541    !!! due to possible change of rot_cmd_max after reading rotation map
542    ALLOCATE(rot_cmd(kjpindex,rot_cmd_max),stat=ier)
543    IF (ier/=0) THEN
544        WRITE(numout,*) "ERROR IN ALLOCATION OF rot_cmd: ",ier
545        STOP 'sechiba_initialize rot_cmd allocation'
546    ENDIF
547    rot_cmd(:,:) = 0
548   
549    !! 1.4 Initialize diffusion coefficients
550    CALL diffuco_initialize (kjit,    kjpindex, index,                  &
551                             rest_id, lalo,     neighbours, resolution, &
552                             rstruct, tq_cdrag, tq_cdrag_pft)
553   
554    !! 1.5 Initialize remaining variables of energy budget
555    CALL enerbil_initialize (kjit,     kjpindex,     index,    rest_id,  &
556                             qair,                             &
557                             temp_sol, temp_sol_pft, temp_sol_new, temp_sol_new_pft, tsol_rad, &
558                             evapot,   evapot_corr,  qsurf,    fluxsens, &
559                             fluxlat,  vevapp )
560   
561   
562    !! 1.7 Initialize remaining hydrological variables
563    IF ( .NOT. hydrol_cwrr ) THEN
564       ! 1.7.1 Initialize remaining hydrological variables from Choisnel module (2 soil layers)
565       
566        CALL hydrolc_initialize( kjit,     kjpindex,     index,          rest_id,   &
567                                  veget,     veget_max,  tot_bare_soil,             &
568                                  rsol,      drysoil_frac, snow,                    &
569                                  snow_age,  snow_nobio,   snow_nobio_age, humrel,  &
570                                  vegstress, qsintveg,     shumdiag,       snowrho, &
571                                  snowtemp,  snowgrain,      snowdz,  &
572                                  snowheat ) 
573     
574       evap_bare_lim(:) = -un
575       k_litt(:) = huit
576       
577       ! No specific calculation for shumdiag_perma. We assume it to shumdiag.
578       shumdiag_perma(:,:)=shumdiag(:,:)
579    ELSE
580       !! 1.7.2 Initialize remaining hydrological variables from CWRR module (11 soil layers)
581          CALL hydrol_initialize (kjit,           kjpindex,  index,         rest_id,          &
582                                  njsc,           soiltile,  veget,         veget_max,        &
583                                  humrel,         vegstress, drysoil_frac,                    &
584                                  shumdiag_perma, qsintveg,                        &
585!                                  evap_bare_lim,  evap_bare_lim_pft, snow,      snow_age,      snow_nobio,       &
586                                  evap_bare_lim,  snow,      snow_age,      snow_nobio,       &
587                                  snow_nobio_age, snowrho,   snowtemp,                        &
588                                  snowgrain,      snowdz,    snowheat,                        &
589                                  totfrac_nobio,  precip_rain, precip_snow, returnflow,       &
590                                  reinfiltration, irrigation,vevapwet,                        &
591                                  transpir,       vevapnu,   vevapsno, vevapflo,              &
592                                  floodout,       runoff,    drainage,                        &
593                                  mc_layh,        mcl_layh,  soilmoist,                       &
594                                  mc_layh_s,      mcl_layh_s,                                 &
595!gmjc
596            tmc_topgrass, humcste_use, altmax, swdown) !!xchen add swdown
597!end gmjc
598    ENDIF
599   
600    !! 1.9 Initialize surface parameters (emissivity, albedo and roughness)
601    CALL condveg_initialize (kjit, kjpindex, index, veget, & 
602            veget_max, frac_nobio, totfrac_nobio, &
603            lalo, neighbours, resolution, contfrac, rest_id, &
604            zlev, drysoil_frac, height, snowdz, snowrho, tot_bare_soil, &
605            snow, snow_age, snow_nobio, snow_nobio_age, &
606            temp_air, pb, u, v, lai, &
607            emis, albedo, z0m, z0h, roughheight, roughheight_pft, &
608            frac_snow_veg, frac_snow_nobio)
609   
610    !! 1.10 Initialization of soil thermodynamics
611    CALL vegetation_nobio_to_bg(veget_max, veget_max_bg)
612    IF (hydrol_cwrr) THEN
613       is_crop_soil(:) = .FALSE.
614       DO jv = 1,nvm
615            mc_layh_pft(:,:,jv) = mc_layh_s(:,:,pref_soil_veg(jv))
616            mcl_layh_pft(:,:,jv) = mcl_layh_s(:,:,pref_soil_veg(jv))
617            soilmoist_pft(:,:,jv) = soilmoist(:, :)
618            IF (ok_LAIdev(jv)) THEN
619                is_crop_soil(pref_soil_veg(jv)) = .TRUE.
620            ENDIF
621       ENDDO
622       IF (printlev >=4) THEN
623           WRITE(numout,*) 'xuhui: initialized'
624           WRITE(numout,*) 'is_crop_soil', is_crop_soil
625           WRITE(numout,*) 'soilmoist(1,1)', soilmoist(1,1)
626           WRITE(numout,*) 'soilmoist_pft(1,1,:)', soilmoist_pft(1,1,:)
627           WRITE(numout,*) 'soilmoist(1,nslm)', soilmoist(1,nslm)
628           WRITE(numout,*) 'soilmoist_pft(1,nslm,:)', soilmoist_pft(1,nslm,:)
629       ENDIF
630       CALL thermosoil_initialize (kjit,  kjpindex, lalo,neighbours,resolution,contfrac, rest_id, veget_max_bg, &
631            shumdiag_perma, snow,               thawed_humidity,  soilc_total, &
632            temp_sol_new,   temp_sol_new_pft, depth_organic_soil, stempdiag, &
633            soilcap,     soilcap_pft,    soilflx,   soilflx_pft,         gtemp, &
634            mc_layh,        mcl_layh,           soilmoist,   mc_layh_pft, mcl_layh_pft, soilmoist_pft, njsc, &
635            frac_snow_veg,  frac_snow_nobio,    totfrac_nobio,     &
636            snowdz,         snowrho,            snowtemp,          &
637            lambda_snow,    cgrnd_snow,         dgrnd_snow,       pb)
638    ELSE
639       CALL thermosoilc_initialize (kjit,  kjpindex, lalo, rest_id, veget_max_bg, &
640            snowdz,        shumdiag_perma, snow,  thawed_humidity,  soilc_total, &
641            temp_sol_new,  depth_organic_soil,    stempdiag, &
642            soilcap,       soilflx,               gtemp ,    & 
643            frac_snow_veg, frac_snow_nobio,       totfrac_nobio,  &
644            snowrho,       snowtemp,          &
645            lambda_snow,    cgrnd_snow,         dgrnd_snow,       pb)
646
647    END IF
648
649    !! 1.12 Initialize river routing
650    IF ( river_routing .AND. nbp_glo .GT. 1) THEN
651       !! 1.12.1 Initialize river routing
652       CALL routing_initialize( kjit,        kjpindex,       index,                 &
653            rest_id,     hist_id,        hist2_id,   lalo,      &
654            neighbours,  resolution,     contfrac,   stempdiag, &
655            returnflow,  reinfiltration, irrigation, riverflow, &
656            coastalflow, flood_frac,     flood_res )
657    ELSE
658       !! 1.12.2 No routing, set variables to zero
659       riverflow(:) = zero
660       coastalflow(:) = zero
661       returnflow(:) = zero
662       reinfiltration(:) = zero
663       irrigation(:) = zero
664       flood_frac(:) = zero
665       flood_res(:) = zero
666    ENDIF
667
668    IF (do_fullirr) THEN ! compulsory irrigation
669        CALL restget_p(rest_id, 'irrigation', nbp_glo, 1, 1, kjit, .TRUE., irrigation, "gather", nbp_glo, index_g)
670          IF (.NOT. (MINVAL(irrigation) < MAXVAL(irrigation) .OR. MAXVAL(irrigation) < val_exp)) THEN
671             irrigation(:) = zero
672          ENDIF
673    ENDIF
674   
675    !! 1.13 Write internal variables to output fields
676    z0m_out(:) = z0m(:)
677    z0h_out(:) = z0h(:)
678    emis_out(:) = emis(:) 
679    qsurf_out(:) = qsurf(:)
680
681    !! 2. Output variables only once
682    zmaxh_glo(:) = zmaxh
683    CALL xios_orchidee_send_field("zmaxh",zmaxh_glo)
684
685    IF (printlev_loc>=3) WRITE(numout,*) 'sechiba_initialize done'
686
687  END SUBROUTINE sechiba_initialize
688
689!! ==============================================================================================================================\n
690!! SUBROUTINE   : sechiba_main
691!!
692!>\BRIEF        Main routine for the sechiba module performing three functions:
693!! calculating temporal evolution of all variables and preparation of output and
694!! restart files (during the last call only)
695!!
696!!\n DESCRIPTION : Main routine for the sechiba module.
697!! One time step evolution consists of:
698!! - call sechiba_var_init to do some initialization,
699!! - call slowproc_main to do some daily calculations
700!! - call diffuco_main for diffusion coefficient calculation,
701!! - call enerbil_main for energy budget calculation,
702!! - call hydrolc_main (or hydrol_main) for hydrologic processes calculation,
703!! - call enerbil_fusion : last part with fusion,
704!! - call condveg_main for surface conditions such as roughness, albedo, and emmisivity,
705!! - call thermosoil_main(for cwrr) or thermosoilc_main(for choisnel) for soil thermodynamic calculation,
706!! - call sechiba_end to swap previous to new fields.
707!!
708!! RECENT CHANGE(S): None
709!!
710!! MAIN OUTPUT VARIABLE(S): Hydrological variables (:: coastalflow and :: riverflow),
711!! components of the energy budget (:: tsol_rad, :: vevapp, :: fluxsens,
712!! :: temp_sol_new and :: fluxlat), surface characteristics (:: z0_out, :: emis_out,
713!! :: tq_cdrag and :: albedo) and land use related CO2 fluxes (:: netco2flux and
714!! :: fco2_lu)           
715!!
716!! REFERENCE(S) :
717!!
718!! FLOWCHART    :
719!! \latexonly
720!! \includegraphics[scale = 0.5]{sechibamainflow.png}
721!! \endlatexonly
722!! \n
723!_ ================================================================================================================================
724
725  SUBROUTINE sechiba_main (kjit, kjpij, kjpindex, index, date0, &
726       & ldrestart_read, ldrestart_write, &
727       & lalo, contfrac, neighbours, resolution,&
728       & zlev, u, v, qair, q2m, t2m, temp_air, epot_air, ccanopy, &
729       & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
730       & precip_rain, precip_snow, lwdown, swnet, swdown, coszang, pb, &
731       & vevapp, fluxsens, fluxlat, coastalflow, riverflow, netco2flux, fco2_lu, &
732       & tsol_rad, temp_sol_new, qsurf_out, albedo, emis_out, z0m_out, z0h_out,&
733       & veget_out, lai_out, height_out, &
734       & rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
735
736!! 0.1 Input variables
737   
738    INTEGER(i_std), INTENT(in)                               :: kjit              !! Time step number (unitless)
739    INTEGER(i_std), INTENT(in)                               :: kjpij             !! Total size of the un-compressed grid
740                                                                                  !! (unitless)
741    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size - terrestrial pixels only
742                                                                                  !! (unitless)
743    INTEGER(i_std),INTENT (in)                               :: rest_id           !! _Restart_ file identifier (unitless)
744    INTEGER(i_std),INTENT (in)                               :: hist_id           !! _History_ file identifier (unitless)
745    INTEGER(i_std),INTENT (in)                               :: hist2_id          !! _History_ file 2 identifier (unitless)
746    INTEGER(i_std),INTENT (in)                               :: rest_id_stom      !! STOMATE's _Restart_ file identifier
747                                                                                  !! (unitless)
748    INTEGER(i_std),INTENT (in)                               :: hist_id_stom      !! STOMATE's _History_ file identifier
749                                                                                  !! (unitless)
750    INTEGER(i_std),INTENT(in)                                :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file
751                                                                                  !! identifier (unitless)
752    REAL(r_std), INTENT (in)                                 :: date0             !! Initial date (??unit??)
753    LOGICAL, INTENT(in)                                      :: ldrestart_read    !! Logical for _restart_ file to read
754                                                                                  !! (true/false)
755    LOGICAL, INTENT(in)                                      :: ldrestart_write   !! Logical for _restart_ file to write
756                                                                                  !! (true/false)
757    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)          :: lalo              !! Geographic coordinates (latitude,longitude)
758                                                                                  !! for grid cells (degrees)
759    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: contfrac          !! Fraction of continent in the grid
760                                                                                  !! (unitless, 0-1)
761    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index             !! Indices of the pixels on the map.
762                                                                                  !! Sechiba uses a reduced grid excluding oceans
763                                                                                  !! ::index contains the indices of the
764                                                                                  !! terrestrial pixels only! (unitless)
765    INTEGER(i_std), DIMENSION(kjpindex,NbNeighb), INTENT(in) :: neighbours        !! Neighboring grid points if land!(unitless)
766    REAL(r_std), DIMENSION (kjpindex,2), INTENT(in)          :: resolution        !! Size in x and y of the grid (m)
767   
768    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u                 !! Lowest level wind speed in direction u
769                                                                                  !! @tex $(m.s^{-1})$ @endtex
770    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: v                 !! Lowest level wind speed in direction v
771                                                                                  !! @tex $(m.s^{-1})$ @endtex
772    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: zlev              !! Height of first layer (m)
773    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair              !! Lowest level specific humidity
774                                                                                  !! @tex $(kg kg^{-1})$ @endtex
775    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q2m               !! 2m specific humidity
776                                                                                  !! @tex $(kg kg^{-1})$ @endtex
777    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: t2m               !! 2m air temperature (K)
778    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: precip_rain       !! Rain precipitation
779                                                                                  !! @tex $(kg m^{-2})$ @endtex
780    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: precip_snow       !! Snow precipitation
781                                                                                  !! @tex $(kg m^{-2})$ @endtex
782    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: lwdown            !! Down-welling long-wave flux
783                                                                                  !! @tex $(W m^{-2})$ @endtex
784    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: coszang           !! Cosine of the solar zenith angle (unitless)
785    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swnet             !! Net surface short-wave flux
786                                                                                  !! @tex $(W m^{-2})$ @endtex
787    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swdown            !! Down-welling surface short-wave flux
788                                                                                  !! @tex $(W m^{-2})$ @endtex
789    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_air          !! Air temperature (K)
790    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: epot_air          !! Air potential energy (??J)
791    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ccanopy           !! CO2 concentration in the canopy (ppm)
792    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petAcoef          !! Coefficients A for T from the Planetary
793                                                                                  !! Boundary Layer
794    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqAcoef          !! Coefficients A for q from the Planetary
795                                                                                  !! Boundary Layer
796    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petBcoef          !! Coefficients B for T from the Planetary
797                                                                                  !! Boundary Layer
798    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqBcoef          !! Coefficients B for q from the Planetary
799                                                                                  !! Boundary Layer
800    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: pb                !! Surface pressure (hPa)
801
802
803!! 0.2 Output variables
804
805    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: coastalflow       !! Outflow on coastal points by small basins.
806                                                                                  !! This is the water which flows in a disperse
807                                                                                  !! way into the ocean
808                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
809    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: riverflow         !! Outflow of the major rivers.
810                                                                                  !! The flux will be located on the continental
811                                                                                  !! grid but this should be a coastal point 
812                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
813    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: tsol_rad          !! Radiative surface temperature
814                                                                                  !! @tex $(W m^{-2})$ @endtex
815    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapp            !! Total of evaporation
816                                                                                  !! @tex $(kg m^{-2} days^{-1})$ @endtex
817   
818    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qsurf_out         !! Surface specific humidity
819                                                                                  !! @tex $(kg kg^{-1})$ @endtex
820    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0m_out           !! Surface roughness momentum (output diagnostic, m)
821    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0h_out           !! Surface roughness heat (output diagnostic, m)
822    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out)         :: albedo            !! Surface albedo for visible and near-infrared
823                                                                                  !! (unitless, 0-1)
824    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxsens          !! Sensible heat flux
825                                                                                  !! @tex $(W m^{-2})$ @endtex
826    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxlat           !! Latent heat flux
827                                                                                  !! @tex $(W m^{-2})$ @endtex
828    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: emis_out          !! Emissivity (unitless)
829    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: netco2flux        !! Sum CO2 flux over PFTs
830                                                                                  !! ??(gC m^{-2} s^{-1})??
831    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fco2_lu           !! Land Cover Change CO2 flux
832                                                                                  !! ??(gC m^{-2} s^{-1})??
833    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: veget_out         !! Fraction of vegetation type (unitless, 0-1)
834    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: lai_out           !! Leaf area index (m^2 m^{-2})
835    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: height_out        !! Vegetation Height (m)
836
837!@chy: these lines are not needed. To be declared as "out", it must be passed as
838!parameters of sechiba_main, but it is not the case here.           
839    !REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: lai_cl1           !! Land Cover Change CO2 flux
840    !REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (put)       :: lai_cl2           !! Land Cover Change CO2 flux
841    !REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: lai_cl3           !! Land Cover Change CO2 flux
842    !REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: lai_cl4           !! Land Cover Change CO2 flux
843    !!REAL(r_std),DIMENSION (kjpindex), INTENT (out)            :: dry_index       !! ratio of
844
845
846!! 0.3 Modified
847
848    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: tq_cdrag          !! Surface drag coefficient  (-)
849    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: temp_sol_new      !! New ground temperature (K)
850
851!! 0.4 local variables
852
853    INTEGER(i_std)                                           :: ji, jv            !! Index (unitless)
854    INTEGER(i_std)                                           :: jsrc,jtar     !! Index (unitless)
855    INTEGER(i_std)                                           :: stveg, edveg
856    REAL(r_std), DIMENSION(kjpindex)                         :: histvar           !! Computations for history files (unitless)
857    REAL(r_std), DIMENSION(kjpindex,nlut)                    :: histvar2          !! Computations for history files (unitless)
858    CHARACTER(LEN=80)                                        :: var_name          !! To store variables names for I/O (unitless)
859    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treefrac      !! Total fraction occupied by trees (0-1, uniless)
860    REAL(r_std), DIMENSION(nvm)                              :: temp_veget_max
861    REAL(r_std), DIMENSION(nvm)                              :: old_veget_max
862    REAL(r_std), DIMENSION(nvm,nvm)                          :: matrix_rot     !! temporary matrix for rotation
863    REAL(r_std)                                              :: rotprc
864    REAL(r_std), DIMENSION(kjpindex)                         :: sum_grassfracC3   !! Total fraction occupied by C3 grasses (0-1, unitless)
865    REAL(r_std), DIMENSION(kjpindex)                         :: sum_grassfracC4   !! Total fraction occupied by C4 grasses (0-1, unitless)
866    REAL(r_std), DIMENSION(kjpindex)                         :: sum_cropfracC3    !! Total fraction occupied by C3 crops (0-1, unitess)
867    REAL(r_std), DIMENSION(kjpindex)                         :: sum_cropfracC4    !! Total fraction occupied by C4 crops (0-1, unitess)
868    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treeFracNdlEvg!! Total fraction occupied by treeFracNdlEvg (0-1, unitess)
869    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treeFracBdlEvg!! Total fraction occupied by treeFracBdlEvg (0-1, unitess)
870    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treeFracNdlDcd!! Total fraction occupied by treeFracNdlDcd (0-1, unitess)
871    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treeFracBdlDcd!! Total fraction occupied by treeFracBdlDcd (0-1, unitess)
872    REAL(r_std), DIMENSION(kjpindex,nsnow)                   :: snowliq           !! Liquid water content (m)
873    REAL(r_std), DIMENSION(kjpindex)                         :: snow_age_diag     !! Only for diag, contains xios_default_val
874    REAL(r_std), DIMENSION(kjpindex,nnobio)                  :: snow_nobio_age_diag !! Only for diag, contains xios_default_val
875    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
876    REAL(r_std), DIMENSION(kjpindex)                         :: grndflux          !! Net energy into soil (W/m2)
877    REAL(r_std),DIMENSION (kjpindex)                         :: wspeed            !! Lowest level Wind speed
878!    REAL(r_std), DIMENSION(kjpindex)                         :: temp_transpot_agr !! potential transpiration by crops need irrigation
879    REAL(r_std)                                              :: temp_irrig_need
880    REAL(r_std)                                              :: tempfrac
881    INTEGER(i_std)                                           :: ier, k
882    REAL(r_std), DIMENSION(kjpindex,nlut)                    :: gpplut            !! GPP on landuse tile, only for diagnostics
883
884    REAL(r_std), DIMENSION(kjpindex,nslm,nvm)                :: mc_layh_pft
885    REAL(r_std), DIMENSION(kjpindex,nslm,nvm)                :: mcl_layh_pft
886    REAL(r_std), DIMENSION(kjpindex,nslm,nvm)                :: soilmoist_pft
887
888    REAL(r_std), DIMENSION(nvm)                              :: temp_temp_sol, temp_soilcap, temp_soilflx
889    REAL(r_std), DIMENSION(nvm)                              :: dilu_temp, dilu_soilcap, dilu_soilflx
890    !!@xchen use var_temp to add FAPAR, PAR_xc, LUE and PC in Wu's paper       
891    REAL(r_std),DIMENSION (kjpindex,nvm)     :: var_temp1                       
892    REAL(r_std),DIMENSION (kjpindex,nvm)     :: var_temp2                       
893    REAL(r_std),DIMENSION (kjpindex,nvm)     :: var_temp3                       
894    REAL(r_std),DIMENSION (kjpindex,nvm)     :: var_temp4 
895    REAL(r_std),DIMENSION (kjpindex,nvm)     :: var_temp5                       
896    REAL(r_std),DIMENSION (kjpindex,nvm)     :: var_temp6 
897    REAL(r_std),DIMENSION (kjpindex,nvm)     :: var_temp10                       
898    REAL(r_std),DIMENSION (kjpindex,nvm)     :: var_temp9                       
899    REAL(r_std),DIMENSION (kjpindex,nvm)     :: var_temp8                       
900    REAL(r_std),DIMENSION (kjpindex,nvm)     :: var_temp7                     
901
902    REAL(r_std), DIMENSION(kjpindex,nvm)                     :: veget_max_bg !! No nobio. All placed in bareground pft
903!_ ================================================================================================================================
904
905    IF (printlev_loc>=3) WRITE(numout,*) 'Start sechiba_main kjpindex =',kjpindex
906
907    !! 1. Initialize variables at each time step
908    CALL sechiba_var_init (kjpindex, rau, pb, temp_air) 
909
910    !! 2. Compute diffusion coefficients
911    CALL diffuco_main (kjit, kjpindex, index, indexveg, indexlai, indexlai0, u, v, &
912         & zlev, z0m, z0h, roughheight, roughheight_pft, temp_sol, temp_sol_pft, temp_air, temp_growth, rau, tq_cdrag, tq_cdrag_pft, &
913         & qsurf, qair, q2m, t2m, pb ,  &
914         & rsol, evap_bare_lim, evapot, evapot_corr, snow, flood_frac, flood_res, frac_nobio, snow_nobio, totfrac_nobio, &
915         & swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai,&
916!!@xchen
917         & laitop, laidow, VPD, &
918         & qsintveg, qsintmax, assim_param, &
919         & vbeta, vbeta_pft, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta4_pft, vbeta5, gsmean, rveget, rstruct, cimean, gpp, &
920         & lalo, neighbours, resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
921         & hist_id, hist2_id)
922   
923    !! 3. Compute energy balance
924    CALL enerbil_main (kjit, kjpindex, &
925         & index, indexveg, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, &
926         & qair, peqAcoef, peqBcoef, pb, rau, vbeta, vbeta_pft, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta4_pft,  vbeta5, &
927         & emis, soilflx, soilflx_pft,  soilcap, soilcap_pft, tq_cdrag, tq_cdrag_pft, veget_max, humrel, fluxsens, fluxlat, &
928         & vevapp, transpir, transpot, vevapnu, vevapnu_pft, vevapwet, vevapsno, vevapflo, temp_sol, temp_sol_pft, tsol_rad, &
929         & temp_sol_new, temp_sol_new_pft, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id,&
930         & precip_rain,snowdz,pgflux,temp_sol_add)
931
932   
933    !! 4. Compute hydrology
934    IF ( .NOT. hydrol_cwrr ) THEN
935       ! 4.1 Water balance from Choisnel module (2 soil layers)
936       vegstress_old = vegstress
937       CALL hydrolc_main (kjit, kjpindex, index, indexveg, &
938            & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, contfrac, vevapwet, veget, veget_max,&
939            & qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age,&
940            & tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, vegstress_old, transpot, humrel, &
941            & vegstress, rsol, drysoil_frac, evapot, evapot_corr, flood_frac, flood_res, shumdiag, litterhumdiag, &
942            & soilcap, rest_id, hist_id, hist2_id, soil_deficit, &
943            & temp_air, pb, u, v, pgflux, &
944            & snowrho,snowtemp,snowgrain,snowdz,snowheat,snowliq, &
945            & grndflux, gtemp, tot_bare_soil, soilflxresid,       &
946              lambda_snow, cgrnd_snow, dgrnd_snow, temp_sol_add)
947
948       evap_bare_lim(:) = -un
949       k_litt(:) = huit
950       
951       ! No specific calculation for shumdiag_perma. We assume it to shumdiag.
952       shumdiag_perma(:,:)=shumdiag(:,:)
953    ELSE
954       !! 4.1 Water balance from CWRR module (11 soil layers)
955       vegstress_old = vegstress
956       CALL hydrol_main (kjit, kjpindex, &
957            & index, indexveg, indexsoil, indexlayer, indexnslm, &
958            & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, &
959            & qsintmax, qsintveg, vevapnu, vevapnu_pft, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age,  &
960            & tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, vegstress_old, transpot, &
961!            & humrel, vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, evap_bare_lim_pft, flood_frac, flood_res, &
962            & humrel, vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, flood_frac, flood_res, &
963            & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope,&
964            & rest_id, hist_id, hist2_id, soil_deficit, is_crop_soil, &
965            & contfrac, stempdiag, &
966            & temp_air, pb, u, v, tq_cdrag, pgflux, &
967            & snowrho, snowtemp, snowgrain, snowdz, snowheat, snowliq, &
968            & grndflux, gtemp, tot_bare_soil, &
969            & soilflxresid, mc_layh, mcl_layh, soilmoist, mc_layh_s, mcl_layh_s, &
970            & lambda_snow, cgrnd_snow, dgrnd_snow, frac_snow_veg, temp_sol_add, &
971!gmjc
972            & tmc_topgrass, humcste_use, swdown) !! xchen add swdown
973!end gmjc
974
975       rsol(:) = -un
976
977    ENDIF
978         
979    !! 5. Compute remaining components of the energy balance
980    IF ( .NOT. ok_explicitsnow ) THEN 
981!!! this means that my code is not compatible with explicit snow due to missing
982!soilcap_pft and temp_sol_new_pft when explicitsnow activated
983       CALL enerbil_fusion (kjpindex, tot_melt, soilcap, soilcap_pft, snowdz, &
984            temp_sol_new, temp_sol_new_pft, fusion)
985    END IF
986
987    !! 6. Compute surface variables (emissivity, albedo and roughness)
988    CALL condveg_main (kjit, kjpindex, index,&
989           lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
990           zlev, snow, snow_age, snow_nobio, snow_nobio_age, tot_bare_soil, &
991           temp_air, pb, u, v, lai, &
992           drysoil_frac, height, snowdz, snowrho, &
993           emis, albedo, frac_snow_veg, frac_snow_nobio, z0m, z0h, roughheight, roughheight_pft, &
994           rest_id, hist_id, hist2_id)
995
996    !! 7. Compute soil thermodynamics
997    CALL vegetation_nobio_to_bg(veget_max, veget_max_bg)
998    IF (hydrol_cwrr) THEN     
999       !!!!! transfer  mc_layh_s,      mcl_layh_s, soilmoist_s
1000       !!!!!       to  mc_layh_pft,    mcl_layh_pft, soilmoist_pft
1001       DO jv = 1,nvm
1002            mc_layh_pft(:,:,jv) = mc_layh_s(:,:,pref_soil_veg(jv))
1003            mcl_layh_pft(:,:,jv) = mcl_layh_s(:,:,pref_soil_veg(jv))
1004            soilmoist_pft(:,:,jv) = soilmoist(:,:)
1005       ENDDO
1006!       WRITE(numout,*) 'xuhui: before thermosoil_main'
1007
1008       CALL thermosoil_main (kjit, kjpindex, &
1009            index, indexgrnd, &
1010            temp_sol_new, temp_sol_new_pft, snow, soilcap, soilcap_pft, soilflx, soilflx_pft, shumdiag_perma, stempdiag, &
1011            ptnlev1, hist_id, hist2_id, &
1012            snowdz, snowrho, snowtemp, gtemp, pb, &
1013            mc_layh, mcl_layh, soilmoist, mc_layh_pft, mcl_layh_pft, soilmoist_pft,  njsc, &
1014            thawed_humidity, depth_organic_soil, heat_Zimov, tdeep, hsdeep,&
1015            soilc_total, veget_max_bg, &
1016            frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
1017            lambda_snow, cgrnd_snow, dgrnd_snow) 
1018    ELSE
1019       CALL thermosoilc_main (kjit, kjpindex, &
1020            index, indexgrnd,indexnslm, &
1021            temp_sol_new, snow, soilcap, soilflx, shumdiag_perma, stempdiag, &
1022            ptnlev1, hist_id, hist2_id, &
1023            snowdz, snowrho, snowtemp, gtemp, pb, &
1024            thawed_humidity, depth_organic_soil, heat_Zimov, tdeep, hsdeep,&
1025            soilc_total, veget_max_bg, &
1026            frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
1027            lambda_snow, cgrnd_snow, dgrnd_snow)
1028    END IF
1029
1030
1031    !! 8. Compute river routing
1032    IF ( river_routing .AND. nbp_glo .GT. 1) THEN
1033       !! 8.1 River routing
1034       CALL routing_main (kjit, kjpindex, index, &
1035            & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget, veget_max, soil_deficit, floodout, runoff, &
1036            & drainage, transpot, evapot_corr, vegstress, precip_rain, humrel, k_litt, flood_frac, flood_res, &
1037            & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
1038    ELSE
1039       !! 8.2 No routing, set variables to zero
1040       riverflow(:) = zero
1041       coastalflow(:) = zero
1042       returnflow(:) = zero
1043       reinfiltration(:) = zero
1044!       irrigation(:) = zero
1045       flood_frac(:) = zero
1046       flood_res(:) = zero
1047
1048       IF ( do_fullirr ) THEN
1049           DO ji = 1, kjpindex
1050               temp_irrig_need = zero
1051               DO jv = 2,nvm
1052                   IF ( veget_max(ji,jv) .GT. 0 ) THEN
1053                       tempfrac = veget(ji,jv)/veget_max(ji,jv)
1054                       IF (tempfrac .LE. 0) THEN
1055                           tempfrac = 0
1056                       ENDIF
1057                       IF ( ok_LAIdev(jv) .AND. (vegstress(ji,jv) .LT. irrig_threshold(jv))  ) THEN
1058                           IF (irrig_drip) THEN
1059                                  temp_irrig_need = temp_irrig_need + irrig_frac(ji) * &
1060                                                    & MIN( irrig_dosmax, irrig_fulfill(jv) * MAX( zero, &
1061    !                                                & transpot(ji,jv)*tempfrac + evapot_corr(ji)*(1-tempfrac) - precip_rain(ji) ) ) &
1062                                                    & transpot(ji,jv)*tempfrac + evapot(ji)*(1-tempfrac) - precip_rain(ji) ) ) &
1063                                                    & * veget_max(ji,jv)
1064                                                    !!!! reconsider if evapot or evapot_corr to be used as irrigation demand
1065                                                    !!!! this also affects the calc of irrig_demand_ratio in hydrol.f90
1066                                                    !!!!! consider adding re-infiltration into this formula, xuhui
1067                           ELSE ! flooding
1068                                 temp_irrig_need = temp_irrig_need + irrig_frac(ji) * &
1069                                                    & MIN( irrig_dosmax, MAX( zero, &
1070                                                    & soil_deficit(ji,jv) )) * veget_max(ji,jv)
1071                                                    !!!!! consider adding re-infiltration into this formula, xuhui
1072                           ENDIF
1073                       ENDIF
1074                   ENDIF
1075               ENDDO
1076               IF (temp_irrig_need .LT. 0) THEN
1077                   !write(*,*) 'sechiba irrigation: negative irrigation need'
1078                   temp_irrig_need  = 0
1079               ENDIF
1080               irrigation(ji) = temp_irrig_need
1081           ENDDO
1082       ELSE
1083            irrigation(:) = zero
1084       ENDIF
1085
1086       CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba)
1087       CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba)
1088    ENDIF
1089
1090    !! 9. Compute slow processes (i.e. 'daily' and annual time step)
1091       !spitfire
1092       ! Compute wind speed for fire model
1093       DO ji = 1, kjpindex
1094          wspeed(ji) = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji)))
1095       ENDDO
1096       !endspit
1097
1098    !! 2.9 Compute slow processes (i.e. 'daily' and annual time step)
1099    ! ::ok_co2 and ::ok_stomate are flags that determine whether the
1100    ! forcing files are written.
1101    CALL slowproc_main (kjit, kjpij, kjpindex, date0, &
1102         index, indexveg, lalo, neighbours, resolution, contfrac, soiltile, fraclut, nwdFraclut, &
1103         temp_air, temp_sol, stempdiag, &
1104         vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, &
1105         !spitfire
1106         wspeed, &
1107         !endspit
1108         gpp, &
1109         deadleaf_cover, &
1110         assim_param, &
1111         lai,&
1112!!@xchen lai_cl1-4 of leaf age class 1-4
1113         lai_cl1, lai_cl2,lai_cl3,lai_cl4,laitop,laidow,VPD,ave_VPD,&
1114         vcmax_cl1, vcmax_cl2,vcmax_cl3,vcmax_cl4,&
1115         frac_age, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
1116         rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
1117         co2_flux, fco2_lu, temp_growth,&
1118         swdown, evapot_corr, & !added for crops, xuhui
1119         tdeep, hsdeep, snow, heat_Zimov, pb, &
1120         sfluxCH4_deep, sfluxCO2_deep, &
1121         thawed_humidity, depth_organic_soil, zz_deep, zz_coef_deep, &
1122         soilc_total,snowdz,snowrho, tot_bare_soil, f_rot_sech, rot_cmd, &
1123!gmjc
1124         tmc_topgrass,humcste_use,altmax)
1125!end gmjc
1126
1127    !! 9.2 Compute global CO2 flux
1128    netco2flux(:) = zero
1129    DO jv = 2,nvm
1130       netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*veget_max(:,jv)
1131    ENDDO
1132    !! 9.3 crop rotation
1133    IF (ok_rotate) THEN
1134       DO ji = 1,kjpindex
1135
1136         IF (f_rot_sech(ji)) THEN
1137            matrix_rot(:,:) = zero
1138            !!! do all the rotation command
1139            k = 1
1140            DO WHILE ( (rot_cmd(ji,k) .GT. 0) .AND. (k .LE. rot_cmd_max) )
1141                CALL sechiba_get_cmd(rot_cmd(ji,k), stveg, edveg, rotprc)
1142                matrix_rot(stveg,edveg) = rotprc
1143                k = k+1
1144            ENDDO
1145            ! The rotation command in rot_cmd has been well checked by
1146            ! in stomate_rotate, no further need to check here, xuhui
1147            !!! converting veget_max
1148            temp_veget_max = veget_max(ji,:)
1149            old_veget_max = veget_max(ji,:)
1150            DO jsrc = 1, nvm
1151                DO jtar = 1,nvm
1152                    IF ( (matrix_rot(jsrc,jtar) .GT. zero) .AND. (matrix_rot(jsrc,jtar) .LT. 1.0+min_sechiba) )  THEN
1153                        temp_veget_max(jsrc) = temp_veget_max(jsrc) - matrix_rot(jsrc,jtar) * veget_max(ji,jsrc)
1154                        temp_veget_max(jtar) = temp_veget_max(jtar) + matrix_rot(jsrc,jtar) * veget_max(ji,jsrc)
1155                    ENDIF
1156                ENDDO
1157            ENDDO
1158            WHERE (temp_veget_max(:) .LT. min_sechiba) 
1159                temp_veget_max(:) = zero
1160            ENDWHERE
1161            veget_max(ji,:) = temp_veget_max
1162           
1163            !!! transfering soil heat and water storage
1164            !!!!!!! veget & soiltile simple update
1165            veget(ji,1) = veget_max(ji,1)
1166            DO jv = 2,nvm
1167                veget(ji,jv) = veget_max(ji,jv) * (un - exp( - lai(ji,jv) * ext_coeff(jv) ) )
1168            ENDDO
1169
1170            soiltile(ji,:) = zero
1171            soiltile(ji,1) = totfrac_nobio(ji)
1172            DO jsrc = 1,nvm
1173                jv = pref_soil_veg(jsrc)
1174                soiltile(ji,jv) = soiltile(ji,jv) + veget_max(ji,jsrc)
1175            ENDDO
1176            IF ((SUM(soiltile(ji,:)) .LT. 1.0-min_sechiba) .OR.  &
1177                (SUM(soiltile(ji,:)) .GT. 1.0+min_sechiba) ) THEN
1178                WRITE(numout,*) 'ji, soiltile(ji,:)', ji, soiltile(ji,:)
1179                STOP 'sechiba_rotation: sum of soiltile not equal to 1'
1180            ENDIF
1181            !!! soil water storage
1182            CALL hydrol_rotation_update (ji, kjpindex, matrix_rot, old_veget_max, veget_max, soiltile, qsintveg)
1183
1184            CALL thermosoil_rotation_update (ji, kjpindex, matrix_rot, old_veget_max) !ptn, cgrnd, dgrnd
1185
1186            !!! sechiba thermol variables: temp_sol_new, soilcap_pft, soilflx_pft,
1187            temp_temp_sol = temp_sol_new_pft(ji,:)
1188            temp_soilcap = soilcap_pft(ji,:)
1189            temp_soilflx = soilflx_pft(ji,:)
1190
1191            DO jtar = 1,nvm
1192                dilu_temp(:) = zero
1193                dilu_soilcap(:) = zero
1194                dilu_soilflx(:) = zero
1195                IF ( SUM(matrix_rot(:,jtar)) .GT. min_sechiba ) THEN
1196                    DO jsrc = 1,nvm
1197                        IF ( matrix_rot(jsrc,jtar) .GT. min_sechiba ) THEN
1198                            dilu_temp(jsrc) = temp_temp_sol(jsrc)
1199                            dilu_soilcap(jsrc) = temp_soilcap(jsrc)
1200                            dilu_soilflx(jsrc) = temp_soilflx(jsrc)
1201                        ENDIF
1202                    ENDDO
1203                    temp_sol_new_pft(ji,jtar) = temp_temp_sol(jtar) * old_veget_max(jtar) * (1.0 - SUM(matrix_rot(jtar,:)))
1204                    soilcap_pft(ji,jtar) = temp_soilcap(jtar) * old_veget_max(jtar) * (1.0 - SUM(matrix_rot(jtar,:)))
1205                    soilflx_pft(ji,jtar) = temp_soilflx(jtar) * old_veget_max(jtar) * (1.0 - SUM(matrix_rot(jtar,:)))
1206                    DO jsrc = 1,nvm
1207                        temp_sol_new_pft(ji,jtar) = temp_sol_new_pft(ji,jtar) + old_veget_max(jsrc) * matrix_rot(jsrc,jtar) * dilu_temp(jsrc)
1208                        soilcap_pft(ji,jtar) = soilcap_pft(ji,jtar) + old_veget_max(jsrc) * matrix_rot(jsrc,jtar) * dilu_soilcap(jsrc)
1209                        soilflx_pft(ji,jtar) = soilflx_pft(ji,jtar) + old_veget_max(jsrc) * matrix_rot(jsrc,jtar) * dilu_soilflx(jsrc)
1210                    ENDDO
1211                    temp_sol_new_pft(ji,jtar) = temp_sol_new_pft(ji,jtar) / veget_max(ji,jtar)
1212                    soilcap_pft(ji,jtar) = soilcap_pft(ji,jtar) / veget_max(ji,jtar)
1213                    soilflx_pft(ji,jtar) = soilflx_pft(ji,jtar) / veget_max(ji,jtar)
1214                ENDIF
1215            ENDDO
1216            !The old temp_sol_new_pft, soilcap_pft, soilflx_pft should be
1217            !maintained as the old value
1218
1219            IF (printlev>=4) THEN
1220                WRITE(numout,*) 'xuhui: debug for sechiba_rotation, ji', ji
1221                WRITE(numout,*) 'rot_cmd(ji,:)', rot_cmd(ji,:)
1222                WRITE(numout,*) 'temp_temp_sol:', temp_temp_sol
1223                WRITE(numout,*) 'temp_sol_new_pft(ji,:)', temp_sol_new_pft(ji,:)
1224            ENDIF
1225
1226            !!! remove the executed rotation commands
1227            f_rot_sech(ji) = .FALSE.
1228            rot_cmd(ji,:) = 0
1229         ENDIF  ! f_rot_sech(ji)
1230       ENDDO ! ji       
1231    ENDIF !ok_rotate
1232    !!!!! end rotation, xuhui
1233 
1234    !! 10. Update the temperature (temp_sol) with newly computed values
1235    CALL sechiba_end (kjpindex, temp_sol_new, temp_sol_new_pft, temp_sol, temp_sol_pft)
1236    !WRITE(numout,*) 'zd sechiba_end 2 ','snowtemp(1,:)',snowtemp(1,:)
1237
1238   
1239    !! 11. Write internal variables to output fields
1240    z0m_out(:) = z0m(:)
1241    z0h_out(:) = z0h(:)
1242    emis_out(:) = emis(:)
1243    qsurf_out(:) = qsurf(:)
1244    veget_out(:,:)  = veget(:,:)
1245    lai_out(:,:)    = lai(:,:)
1246    height_out(:,:) = height(:,:)
1247 
1248    !! 12. Write global variables to history files
1249    sum_treefrac(:) = zero
1250    sum_grassfracC3(:) = zero
1251    sum_grassfracC4(:) = zero
1252    sum_cropfracC3(:) = zero
1253    sum_cropfracC4(:) = zero
1254    sum_treeFracNdlEvg(:) = zero
1255    sum_treeFracBdlEvg(:) = zero
1256    sum_treeFracNdlDcd(:) = zero
1257    sum_treeFracBdlDcd(:) = zero
1258
1259    DO jv = 2, nvm 
1260       IF (is_tree(jv) .AND. natural(jv)) THEN
1261          sum_treefrac(:) = sum_treefrac(:) + veget_max(:,jv)
1262       ELSE IF ((.NOT. is_tree(jv))  .AND. natural(jv)) THEN
1263          ! Grass
1264          IF (is_c4(jv)) THEN
1265             sum_grassfracC4(:) = sum_grassfracC4(:) + veget_max(:,jv)
1266          ELSE
1267             sum_grassfracC3(:) = sum_grassfracC3(:) + veget_max(:,jv)
1268          END IF
1269       ELSE 
1270          ! Crop and trees not natural
1271          IF (is_c4(jv)) THEN
1272             sum_cropfracC4(:) = sum_cropfracC4(:) + veget_max(:,jv)
1273          ELSE
1274             sum_cropfracC3(:) = sum_cropfracC3(:) + veget_max(:,jv)
1275          END IF
1276       ENDIF
1277
1278       IF (is_tree(jv)) THEN
1279          IF (is_evergreen(jv)) THEN
1280             IF (is_needleleaf(jv)) THEN
1281                ! Fraction for needleleaf evergreen trees (treeFracNdlEvg)
1282                sum_treeFracNdlEvg(:) = sum_treeFracNdlEvg(:) + veget_max(:,jv)
1283             ELSE
1284                ! Fraction for broadleaf evergreen trees (treeFracBdlEvg)
1285                sum_treeFracBdlEvg(:) = sum_treeFracBdlEvg(:) + veget_max(:,jv)
1286             END IF
1287          ELSE IF (is_deciduous(jv)) THEN
1288             IF (is_needleleaf(jv)) THEN
1289                ! Fraction for needleleaf deciduous trees (treeFracNdlDcd)
1290                sum_treeFracNdlDcd(:) = sum_treeFracNdlDcd(:) + veget_max(:,jv)
1291             ELSE 
1292                ! Fraction for broadleafs deciduous trees (treeFracBdlDcd)
1293                sum_treeFracBdlDcd(:) = sum_treeFracBdlDcd(:) + veget_max(:,jv)
1294             END IF
1295          END IF
1296       END IF
1297    ENDDO         
1298
1299    histvar(:)=zero
1300    DO jv = 2, nvm
1301       IF (is_deciduous(jv)) THEN
1302          histvar(:) = histvar(:) + veget_max(:,jv)*100*contfrac
1303       ENDIF
1304    ENDDO
1305    CALL xios_orchidee_send_field("treeFracPrimDec",histvar)
1306
1307    histvar(:)=zero
1308    DO jv = 2, nvm
1309       IF (is_evergreen(jv)) THEN
1310          histvar(:) = histvar(:) + veget_max(:,jv)*100*contfrac
1311       ENDIF
1312    ENDDO
1313    CALL xios_orchidee_send_field("treeFracPrimEver",histvar)
1314
1315    histvar(:)=zero
1316    DO jv = 2, nvm
1317       IF ( .NOT.(is_c4(jv)) ) THEN
1318          histvar(:) = histvar(:) + veget_max(:,jv)*100*contfrac
1319       ENDIF
1320    ENDDO
1321    CALL xios_orchidee_send_field("c3PftFrac",histvar)
1322
1323    histvar(:)=zero
1324    DO jv = 2, nvm
1325       IF ( is_c4(jv) ) THEN
1326          histvar(:) = histvar(:) + veget_max(:,jv)*100*contfrac
1327       ENDIF
1328    ENDDO
1329    CALL xios_orchidee_send_field("c4PftFrac",histvar)
1330
1331    CALL xios_orchidee_send_field("temp_sol_new",temp_sol_new)
1332    CALL xios_orchidee_send_field("fluxsens",fluxsens)
1333    CALL xios_orchidee_send_field("fluxlat",fluxlat)
1334
1335
1336    ! Add XIOS default value where no snow
1337    DO ji=1,kjpindex 
1338       IF (snow(ji) .GT. zero) THEN
1339          snow_age_diag(ji) = snow_age(ji)
1340          snow_nobio_age_diag(ji,:) = snow_nobio_age(ji,:)
1341       
1342          snowage_glob(ji) = snow_age(ji)*frac_snow_veg(ji)*(1-totfrac_nobio(ji)) + &
1343               SUM(snow_nobio_age(ji,:)*frac_snow_nobio(ji,:)*frac_nobio(ji,:))
1344          IF (snowage_glob(ji) .NE. 0) snowage_glob(ji) = snowage_glob(ji) / &
1345               (frac_snow_veg(ji)*(1-totfrac_nobio(ji)) + SUM(frac_snow_nobio(ji,:)*frac_nobio(ji,:)))
1346       ELSE
1347          snow_age_diag(ji) = xios_default_val
1348          snow_nobio_age_diag(ji,:) = xios_default_val
1349          snowage_glob(ji) = xios_default_val
1350       END IF
1351    END DO
1352   
1353    CALL xios_orchidee_send_field("snow",snow)
1354    CALL xios_orchidee_send_field("snowage",snow_age_diag)
1355    CALL xios_orchidee_send_field("snownobio",snow_nobio)
1356    CALL xios_orchidee_send_field("snownobioage",snow_nobio_age_diag)
1357    CALL xios_orchidee_send_field("snowage_glob",snowage_glob)
1358
1359    CALL xios_orchidee_send_field("frac_snow", SUM(frac_snow_nobio,2)*totfrac_nobio+frac_snow_veg*(1-totfrac_nobio))
1360    CALL xios_orchidee_send_field("frac_snow_veg", frac_snow_veg)
1361    CALL xios_orchidee_send_field("frac_snow_nobio", frac_snow_nobio)
1362    CALL xios_orchidee_send_field("reinf_slope",reinf_slope)
1363    CALL xios_orchidee_send_field("njsc",REAL(njsc, r_std))
1364    CALL xios_orchidee_send_field("vegetfrac",veget)
1365    CALL xios_orchidee_send_field("maxvegetfrac",veget_max)
1366    CALL xios_orchidee_send_field("nobiofrac",frac_nobio)
1367    CALL xios_orchidee_send_field("soiltile",soiltile)
1368    CALL xios_orchidee_send_field("rstruct",rstruct)
1369    IF (ok_co2) CALL xios_orchidee_send_field("gpp",gpp/dt_sechiba)
1370
1371    histvar(:)=zero
1372    DO jv = 2, nvm
1373       IF ( .NOT. is_tree(jv) .AND. natural(jv) ) THEN
1374          histvar(:) = histvar(:) + gpp(:,jv)
1375       ENDIF
1376    ENDDO
1377    CALL xios_orchidee_send_field("gppgrass",histvar/dt_sechiba)
1378
1379    histvar(:)=zero
1380    DO jv = 2, nvm
1381       IF ( (.NOT. is_tree(jv)) .AND. (.NOT. natural(jv)) ) THEN
1382          histvar(:) = histvar(:) + gpp(:,jv)
1383       ENDIF
1384    ENDDO
1385    CALL xios_orchidee_send_field("gppcrop",histvar/dt_sechiba)
1386
1387    histvar(:)=zero
1388    DO jv = 2, nvm
1389       IF ( is_tree(jv) ) THEN
1390          histvar(:) = histvar(:) + gpp(:,jv)
1391       ENDIF
1392    ENDDO
1393    CALL xios_orchidee_send_field("gpptree",histvar/dt_sechiba)
1394
1395    CALL xios_orchidee_send_field("nee",co2_flux/dt_sechiba)
1396    CALL xios_orchidee_send_field("drysoil_frac",drysoil_frac)
1397    CALL xios_orchidee_send_field("vevapflo",vevapflo/dt_sechiba)
1398    CALL xios_orchidee_send_field("k_litt",k_litt)
1399    CALL xios_orchidee_send_field("beta",vbeta)
1400    CALL xios_orchidee_send_field("vbeta1",vbeta1)
1401    CALL xios_orchidee_send_field("vbeta2",vbeta2)
1402    CALL xios_orchidee_send_field("vbeta3",vbeta3)
1403    CALL xios_orchidee_send_field("vbeta4",vbeta4)
1404    CALL xios_orchidee_send_field("vbeta5",vbeta5)
1405    CALL xios_orchidee_send_field("gsmean",gsmean)
1406    CALL xios_orchidee_send_field("cimean",cimean)
1407    CALL xios_orchidee_send_field("rveget",rveget)
1408    CALL xios_orchidee_send_field("rsol",rsol)
1409 
1410    histvar(:)=SUM(vevapwet(:,:),dim=2)
1411    CALL xios_orchidee_send_field("evspsblveg",histvar/dt_sechiba)
1412    histvar(:)= vevapnu(:)+vevapsno(:)
1413    CALL xios_orchidee_send_field("evspsblsoi",histvar/dt_sechiba)
1414    histvar(:)=SUM(transpir(:,:),dim=2)
1415    CALL xios_orchidee_send_field("tran",histvar/dt_sechiba)
1416
1417    ! For CMIP6 data request: the following fractions are fractions of the total grid-cell,
1418    ! which explains the multiplication by contfrac
1419    CALL xios_orchidee_send_field("treeFrac",sum_treefrac(:)*100*contfrac(:))
1420    CALL xios_orchidee_send_field("grassFracC3",sum_grassfracC3(:)*100*contfrac(:))
1421    CALL xios_orchidee_send_field("grassFracC4",sum_grassfracC4(:)*100*contfrac(:))
1422    CALL xios_orchidee_send_field("cropFracC3",sum_cropfracC3(:)*100*contfrac(:))
1423    CALL xios_orchidee_send_field("cropFracC4",sum_cropfracC4(:)*100*contfrac(:))
1424    CALL xios_orchidee_send_field("treeFracNdlEvg",sum_treeFracNdlEvg(:)*100*contfrac(:))
1425    CALL xios_orchidee_send_field("treeFracBdlEvg",sum_treeFracBdlEvg(:)*100*contfrac(:))
1426    CALL xios_orchidee_send_field("treeFracNdlDcd",sum_treeFracNdlDcd(:)*100*contfrac(:))
1427    CALL xios_orchidee_send_field("treeFracBdlDcd",sum_treeFracBdlDcd(:)*100*contfrac(:))
1428
1429    histvar(:)=veget_max(:,1)*100*contfrac(:)
1430    CALL xios_orchidee_send_field("baresoilFrac",histvar)
1431    histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)*100*contfrac(:)
1432    CALL xios_orchidee_send_field("residualFrac",histvar)
1433
1434    ! For CMIP6 data request: cnc = canopy cover fraction over land area
1435    histvar(:)=zero
1436    DO jv=2,nvm
1437       histvar(:) = histvar(:) + veget_max(:,jv)*100
1438    END DO
1439    CALL xios_orchidee_send_field("cnc",histvar)
1440   
1441    CALL xios_orchidee_send_field("tsol_rad",tsol_rad-273.15)
1442    CALL xios_orchidee_send_field("qsurf",qsurf)
1443    CALL xios_orchidee_send_field("emis",emis)
1444    CALL xios_orchidee_send_field("z0m",z0m)
1445    CALL xios_orchidee_send_field("z0h",z0h)
1446    CALL xios_orchidee_send_field("roughheight",roughheight)
1447    CALL xios_orchidee_send_field("roughheight_pft",roughheight_pft)
1448    CALL xios_orchidee_send_field("lai",lai)
1449    !!@xchen add PC_wu
1450    !!@xchen add FAPAR=(1-p1)-(1-p2)*exp(-0.5*lai/coszang); p1=0.05 and p2=0 for tower site Santarem km67
1451    !!umol photons/m^2/si
1452    !!lue:0.0088176, 0.02162689,0.034436,0.0211652 for lai_cl1-4                                                         
1453    var_temp1(:,:) = zero                                                       
1454    var_temp2(:,:) = zero                                                       
1455    var_temp3(:,:) = zero                                                       
1456    var_temp4(:,:) = zero
1457    var_temp5(:,:) = zero                                                       
1458    var_temp6(:,:) = zero
1459    var_temp7(:,:) = zero                                                       
1460    var_temp8(:,:) = zero                                                       
1461    var_temp9(:,:) = zero                                                       
1462    var_temp10(:,:) = zero                                                       
1463                                                       
1464                                                                               
1465    DO jv = 1,nvm                                                               
1466       DO ji = 1, kjpindex
1467          IF ( coszang(ji) .NE. 0 ) THEN
1468             var_temp1(ji,jv) = 1-0.05-exp(-0.5*lai(ji,jv)/coszang(ji))                 
1469          ENDIF
1470       ENDDO
1471       !!var_temp2(:,jv) = 0.0088176*lai_cl1(:,jv)+0.02162689*lai_cl2(:,jv)+0.034436*lai_cl3(:,jv)+0.0211652*lai_cl4(:,jv)
1472       var_temp2(:,jv) = 0.0288176*lai_cl1(:,jv)+0.034436*lai_cl2(:,jv)+0.0211652*lai_cl3(:,jv)+0.00211652*lai_cl4(:,jv) 
1473       var_temp5(:,jv) = vcmax_cl1(:,jv)/1210.432*lai_cl1(:,jv)+vcmax_cl2(:,jv)/1146.755*lai_cl2(:,jv)+ &
1474                         & vcmax_cl3(:,jv)/1083.079*lai_cl3(:,jv)+vcmax_cl4(:,jv)/1102.139*lai_cl4(:,jv)
1475 
1476       !!WHERE (var_temp2(:,jv) .GT. zero .AND. var_temp1(:,jv) .GT. zero)     
1477       var_temp3(:,jv) = lai_cl1(:,jv)+lai_cl2(:,jv)+lai_cl3(:,jv)+lai_cl4(:,jv)                       
1478       DO ji = 1, kjpindex
1479          IF (var_temp3(ji,jv) .NE. 0) THEN
1480               var_temp4(ji,jv) =var_temp1(ji,jv)*var_temp2(ji,jv)/var_temp3(ji,jv)
1481               var_temp6(ji,jv) =var_temp1(ji,jv)*var_temp5(ji,jv)/var_temp3(ji,jv)
1482          ENDIF
1483       ENDDO       
1484       var_temp7(:,jv) = vcmax_cl1(:,jv)/1210.432*swdown(:)*W_to_mol*RG_to_PAR*var_temp1(:,jv)
1485       var_temp8(:,jv) = vcmax_cl2(:,jv)/1146.755*swdown(:)*W_to_mol*RG_to_PAR*var_temp1(:,jv) 
1486       var_temp9(:,jv) = vcmax_cl3(:,jv)/1083.079*swdown(:)*W_to_mol*RG_to_PAR*var_temp1(:,jv)
1487       var_temp10(:,jv) = vcmax_cl4(:,jv)/1102.139*swdown(:)*W_to_mol*RG_to_PAR*var_temp1(:,jv) 
1488       !!var_temp4(:,jv) = var_temp1(:,jv)*var_temp2(:,jv)                                             
1489       !!var_temp3(:,jv) = var_temp2(:,jv)*var_temp1(:,jv)                     
1490       !!ENDWHERE                                                               
1491    ENDDO   
1492    CALL xios_orchidee_send_field("PC_wu",var_temp4)
1493    CALL xios_orchidee_send_field("PC_xc",2.5*var_temp6)                           
1494    CALL xios_orchidee_send_field("GPP_CL1",var_temp7) 
1495    CALL xios_orchidee_send_field("GPP_CL2",var_temp8)
1496    CALL xios_orchidee_send_field("GPP_CL3",var_temp9)
1497    CALL xios_orchidee_send_field("GPP_CL4",var_temp10)
1498    CALL xios_orchidee_send_field("GPP_xc",var_temp7+var_temp8+var_temp9+var_temp10)     
1499    CALL histwrite_p(hist_id, 'gpp_cl1', kjit, var_temp7, kjpindex*nvm, indexveg)
1500    CALL histwrite_p(hist_id, 'gpp_cl2', kjit, var_temp8, kjpindex*nvm, indexveg)
1501    CALL histwrite_p(hist_id, 'gpp_cl3', kjit, var_temp9, kjpindex*nvm, indexveg)
1502    CALL histwrite_p(hist_id, 'gpp_cl4', kjit, var_temp10, kjpindex*nvm, indexveg)
1503    CALL histwrite_p(hist_id, 'gpp_xc', kjit, var_temp7+var_temp8+var_temp9+var_temp10, kjpindex*nvm, indexveg)
1504
1505    histvar(:)=zero   
1506    DO ji = 1, kjpindex
1507       IF (SUM(veget_max(ji,:)) > zero) THEN
1508         DO jv=2,nvm
1509            histvar(ji) = histvar(ji) + veget_max(ji,jv)*lai(ji,jv)/SUM(veget_max(ji,:))
1510         END DO
1511       END IF
1512    END DO
1513    !!xchen add P-ET
1514    histvar(:)=zero                                                             
1515    DO jv=1,nvm                                                                 
1516      histvar(:) = histvar(:) + transpir(:,jv)                                 
1517    ENDDO 
1518    !!CALL xios_orchidee_send_field("water_profit",precip_rain-histvar*one_day/dt_sechiba-vevapp*one_day/dt_sechiba)
1519    !!CALL xios_orchidee_send_field("MCWD",evapot_corr-histvar*one_day/dt_sechiba-vevapp*one_day/dt_sechiba)
1520    CALL xios_orchidee_send_field("PET",evapot_corr/dt_sechiba)
1521    CALL xios_orchidee_send_field("LAImean",histvar)
1522    CALL xios_orchidee_send_field("vevapsno",vevapsno/dt_sechiba)
1523    CALL xios_orchidee_send_field("vevapp",vevapp/dt_sechiba)
1524    CALL xios_orchidee_send_field("vevapnu_pft",vevapnu_pft/dt_sechiba)
1525    CALL xios_orchidee_send_field("vevapnu",vevapnu/dt_sechiba)
1526    CALL xios_orchidee_send_field("transpir",transpir*one_day/dt_sechiba)
1527    CALL xios_orchidee_send_field("inter",vevapwet*one_day/dt_sechiba)
1528    CALL xios_orchidee_send_field("Qf",fusion)
1529    IF (.NOT. ( river_routing .AND. nbp_glo .GT. 1) ) THEN
1530        CALL xios_orchidee_send_field("irrigation",irrigation*one_day/dt_sechiba)
1531    ENDIF
1532    histvar(:)=zero
1533    DO jv=1,nvm
1534      histvar(:) = histvar(:) + vevapwet(:,jv)
1535    ENDDO
1536    CALL xios_orchidee_send_field("ECanop",histvar/dt_sechiba)
1537    histvar(:)=zero
1538    DO jv=1,nvm
1539      histvar(:) = histvar(:) + transpir(:,jv)
1540    ENDDO
1541    CALL xios_orchidee_send_field("TVeg",histvar/dt_sechiba)
1542    CALL xios_orchidee_send_field("ACond",tq_cdrag)
1543    CALL xios_orchidee_send_field("ACond_pft",tq_cdrag_pft)
1544
1545
1546    !! Calculate diagnostic variables on Landuse tiles for LUMIP/CMIP6
1547
1548    ! Calculate fraction of landuse tiles related to the whole grid cell
1549    DO jv=1,nlut
1550       histvar2(:,jv) = fraclut(:,jv) * contfrac(:)
1551    END DO
1552    CALL xios_orchidee_send_field("fraclut",histvar2)
1553
1554    CALL xios_orchidee_send_field("nwdFraclut",nwdFraclut(:,:))
1555   
1556    ! Calculate GPP on landuse tiles
1557    ! val_exp is used as missing value where the values are not known i.e. where the tile is not represented
1558    ! or for pasture (id_pst) or urban land (id_urb).
1559    gpplut(:,:)=0
1560    DO jv=1,nvm
1561       IF (natural(jv)) THEN
1562          gpplut(:,id_psl) = gpplut(:,id_psl) + gpp(:,jv)
1563       ELSE
1564          gpplut(:,id_crp) = gpplut(:,id_crp) + gpp(:,jv)
1565       ENDIF
1566    END DO
1567
1568    ! Transform from gC/m2/s into kgC/m2/s
1569    WHERE (fraclut(:,id_psl)>min_sechiba)
1570       gpplut(:,id_psl) = gpplut(:,id_psl)/fraclut(:,id_psl)/1000
1571    ELSEWHERE
1572       gpplut(:,id_psl) = xios_default_val
1573    END WHERE
1574    WHERE (fraclut(:,id_crp)>min_sechiba)
1575       gpplut(:,id_crp) = gpplut(:,id_crp)/fraclut(:,id_crp)/1000
1576    ELSEWHERE
1577       gpplut(:,id_crp) = xios_default_val
1578    END WHERE
1579    gpplut(:,id_pst) = xios_default_val
1580    gpplut(:,id_urb) = xios_default_val
1581
1582    CALL xios_orchidee_send_field("gpplut",gpplut)
1583
1584
1585    IF ( .NOT. almaoutput ) THEN
1586       ! Write history file in IPSL-format
1587       CALL histwrite_p(hist_id, 'beta', kjit, vbeta, kjpindex, index)
1588       CALL histwrite_p(hist_id, 'beta_pft', kjit, vbeta_pft, kjpindex*nvm, indexveg)
1589       CALL histwrite_p(hist_id, 'z0m', kjit, z0m, kjpindex, index)
1590       CALL histwrite_p(hist_id, 'z0h', kjit, z0h, kjpindex, index)
1591       CALL histwrite_p(hist_id, 'soilflx',kjit, soilflx, kjpindex, index)
1592       CALL histwrite_p(hist_id, 'soilcap',kjit, soilcap, kjpindex, index)
1593       CALL histwrite_p(hist_id, 'soilflx_pft',kjit, soilflx_pft, kjpindex, indexveg)
1594       CALL histwrite_p(hist_id, 'soilcap_pft',kjit, soilcap_pft, kjpindex, indexveg)
1595       CALL histwrite_p(hist_id, 'roughheight', kjit, roughheight, kjpindex, index)
1596       CALL histwrite_p(hist_id, 'roughheight_pft', kjit, roughheight_pft, kjpindex, indexveg)
1597       CALL histwrite_p(hist_id, 'temp_sol_pft', kjit, temp_sol_pft, kjpindex, indexveg)
1598       CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1599       CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1600       CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1601       CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1602       !CALL histwrite_p(hist_id, 'lai_cl1', kjit, lai_cl1, kjpindex*nvm, indexveg)
1603       !CALL histwrite_p(hist_id, 'lai_cl2', kjit, lai_cl2, kjpindex*nvm, indexveg)
1604       !CALL histwrite_p(hist_id, 'lai_cl3', kjit, lai_cl3, kjpindex*nvm, indexveg)
1605       !CALL histwrite_p(hist_id, 'lai_cl4', kjit, lai_cl4, kjpindex*nvm, indexveg)
1606       CALL histwrite_p(hist_id, 'subli', kjit, vevapsno, kjpindex, index)
1607       CALL histwrite_p(hist_id, 'evapnu', kjit, vevapnu, kjpindex, index)
1608       CALL histwrite_p(hist_id, 'evapnu_pft', kjit, vevapnu_pft, kjpindex*nvm, indexveg)
1609       CALL histwrite_p(hist_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
1610       CALL histwrite_p(hist_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
1611       CALL histwrite_p(hist_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
1612       CALL histwrite_p(hist_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
1613       CALL histwrite_p(hist_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
1614       CALL histwrite_p(hist_id, 'vbeta4', kjit, vbeta4, kjpindex, index)   
1615       CALL histwrite_p(hist_id, 'vbeta4_pft', kjit, vbeta4_pft, kjpindex*nvm, indexveg)   
1616       CALL histwrite_p(hist_id, 'vbeta5', kjit, vbeta5, kjpindex, index)   
1617       CALL histwrite_p(hist_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
1618       CALL histwrite_p(hist_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
1619       CALL histwrite_p(hist_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
1620
1621       IF ( .NOT. hydrol_cwrr ) THEN
1622          CALL histwrite_p(hist_id, 'rsol', kjit, rsol, kjpindex, index)
1623       ENDIF
1624       CALL histwrite_p(hist_id, 'snow', kjit, snow, kjpindex, index)
1625       CALL histwrite_p(hist_id, 'snowage', kjit, snow_age, kjpindex, index)
1626       CALL histwrite_p(hist_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
1627       CALL histwrite_p(hist_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
1628       IF (.NOT. ( river_routing .AND. nbp_glo .GT. 1) ) THEN 
1629       ! if routing not activated output here
1630       ! otherwise output is inside routing
1631           CALL histwrite_p(hist_id, 'irrigation', kjit, irrigation, kjpindex, index)
1632       ENDIF
1633
1634       IF (ok_explicitsnow) THEN
1635          CALL histwrite_p(hist_id, 'grndflux', kjit, grndflux, kjpindex,index)
1636          CALL histwrite_p(hist_id, 'snowtemp',kjit,snowtemp,kjpindex*nsnow,indexsnow)
1637          CALL histwrite_p(hist_id, 'snowliq', kjit,snowliq,kjpindex*nsnow,indexsnow)
1638          CALL histwrite_p(hist_id, 'snowdz', kjit,snowdz,kjpindex*nsnow,indexsnow)
1639          CALL histwrite_p(hist_id, 'snowrho', kjit,snowrho,kjpindex*nsnow,indexsnow)
1640          CALL histwrite_p(hist_id, 'snowgrain',kjit,snowgrain,kjpindex*nsnow,indexsnow)
1641          CALL histwrite_p(hist_id, 'snowheat',kjit,snowheat,kjpindex*nsnow,indexsnow)
1642       END IF
1643
1644       CALL histwrite_p(hist_id, 'pgflux',kjit,pgflux,kjpindex,index)
1645       CALL histwrite_p(hist_id, 'soiltile',  kjit, soiltile, kjpindex*nstm, indexsoil)
1646       !
1647       IF ( hydrol_cwrr ) THEN
1648          CALL histwrite_p(hist_id, 'soilindex',  kjit, REAL(njsc, r_std), kjpindex, index)
1649          CALL histwrite_p(hist_id, 'reinf_slope',  kjit, reinf_slope, kjpindex, index)
1650          CALL histwrite_p(hist_id, 'k_litt', kjit, k_litt, kjpindex, index)
1651       ENDIF
1652       IF ( river_routing .AND. do_floodplains ) THEN
1653          CALL histwrite_p(hist_id, 'evapflo', kjit, vevapflo, kjpindex, index)
1654          CALL histwrite_p(hist_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
1655       ENDIF
1656       IF ( ok_co2 ) THEN
1657          CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)   
1658          CALL histwrite_p(hist_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
1659          CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)   
1660       ENDIF
1661       IF ( ok_stomate ) THEN
1662          CALL histwrite_p(hist_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)   
1663       ENDIF
1664
1665       histvar(:)=SUM(vevapwet(:,:),dim=2)
1666       CALL histwrite_p(hist_id, 'evspsblveg', kjit, histvar, kjpindex, index)
1667
1668       histvar(:)= vevapnu(:)+vevapsno(:)
1669       CALL histwrite_p(hist_id, 'evspsblsoi', kjit, histvar, kjpindex, index)
1670
1671       histvar(:)=SUM(transpir(:,:),dim=2)
1672       CALL histwrite_p(hist_id, 'tran', kjit, histvar, kjpindex, index)
1673
1674       histvar(:)= sum_treefrac(:)*100*contfrac(:)
1675       CALL histwrite_p(hist_id, 'treeFrac', kjit, histvar, kjpindex, index) 
1676
1677       histvar(:)= (sum_grassfracC3(:)+sum_grassfracC4(:))*100*contfrac(:)
1678       CALL histwrite_p(hist_id, 'grassFrac', kjit, histvar, kjpindex, index) 
1679
1680       histvar(:)= (sum_cropfracC3(:)+sum_cropfracC4(:))*100*contfrac(:)
1681       CALL histwrite_p(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)
1682
1683       histvar(:)=veget_max(:,1)*100*contfrac(:)
1684       CALL histwrite_p(hist_id, 'baresoilFrac', kjit, histvar, kjpindex, index)
1685
1686       histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)*100*contfrac(:)
1687       CALL histwrite_p(hist_id, 'residualFrac', kjit, histvar, kjpindex, index)
1688    ELSE
1689       ! Write history file in ALMA format
1690       CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1691       CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1692       CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1693       CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1694       IF (.NOT. ok_explicitsnow) CALL histwrite_p(hist_id, 'Qf', kjit, fusion, kjpindex, index)
1695       CALL histwrite_p(hist_id, 'ESoil', kjit, vevapnu, kjpindex, index)
1696       CALL histwrite_p(hist_id, 'EWater', kjit, vevapflo, kjpindex, index)
1697       CALL histwrite_p(hist_id, 'SWE', kjit, snow, kjpindex, index)
1698       histvar(:)=zero
1699       DO jv=1,nvm
1700          histvar(:) = histvar(:) + transpir(:,jv)
1701       ENDDO
1702       CALL histwrite_p(hist_id, 'TVeg', kjit, histvar, kjpindex, index)
1703       histvar(:)=zero
1704       DO jv=1,nvm
1705          histvar(:) = histvar(:) + vevapwet(:,jv)
1706       ENDDO
1707       CALL histwrite_p(hist_id, 'ECanop', kjit, histvar, kjpindex, index)
1708       CALL histwrite_p(hist_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
1709       CALL histwrite_p(hist_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
1710       !
1711       CALL histwrite_p(hist_id, 'Z0m', kjit, z0m, kjpindex, index)
1712       CALL histwrite_p(hist_id, 'Z0h', kjit, z0h, kjpindex, index)
1713       CALL histwrite_p(hist_id, 'EffectHeight', kjit, roughheight, kjpindex, index)
1714       !
1715       IF ( river_routing .AND. do_floodplains ) THEN
1716          CALL histwrite_p(hist_id, 'Qflood', kjit, vevapflo, kjpindex, index)
1717          CALL histwrite_p(hist_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
1718       ENDIF
1719       !
1720       IF ( ok_co2 ) THEN
1721          CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)   
1722          CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)   
1723          CALL histwrite_p(hist_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
1724       ENDIF
1725       IF ( ok_stomate ) THEN
1726             CALL histwrite_p(hist_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)   
1727       ENDIF
1728    ENDIF ! almaoutput
1729   
1730    !! 13. Write additional output file with higher frequency
1731    IF ( hist2_id > 0 ) THEN
1732       IF ( .NOT. almaoutput ) THEN
1733          ! Write history file in IPSL-format
1734          CALL histwrite_p(hist2_id, 'tsol_rad', kjit, tsol_rad, kjpindex, index)
1735          CALL histwrite_p(hist2_id, 'qsurf', kjit, qsurf, kjpindex, index)
1736          CALL histwrite_p(hist2_id, 'albedo', kjit, albedo, kjpindex*2, indexalb)
1737          CALL histwrite_p(hist2_id, 'emis', kjit, emis, kjpindex, index)
1738          CALL histwrite_p(hist2_id, 'beta', kjit, vbeta, kjpindex, index)
1739          CALL histwrite_p(hist2_id, 'z0m', kjit, z0m, kjpindex, index)
1740          CALL histwrite_p(hist2_id, 'z0h', kjit, z0h, kjpindex, index)
1741          CALL histwrite_p(hist2_id, 'roughheight', kjit, roughheight, kjpindex, index)
1742          CALL histwrite_p(hist2_id, 'roughheight_pft', kjit, roughheight_pft, kjpindex*nvm, indexveg)
1743          CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1744          CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1745          CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1746          CALL histwrite_p(hist2_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
1747          CALL histwrite_p(hist2_id, 'subli', kjit, vevapsno, kjpindex, index)
1748          IF ( river_routing .AND. do_floodplains ) THEN
1749             CALL histwrite_p(hist2_id, 'vevapflo', kjit, vevapflo, kjpindex, index)
1750             CALL histwrite_p(hist2_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
1751          ENDIF
1752          CALL histwrite_p(hist2_id, 'vevapnu', kjit, vevapnu, kjpindex, index)
1753          CALL histwrite_p(hist2_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
1754          CALL histwrite_p(hist2_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
1755          CALL histwrite_p(hist2_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
1756          CALL histwrite_p(hist2_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
1757          CALL histwrite_p(hist2_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
1758          CALL histwrite_p(hist2_id, 'vbeta4', kjit, vbeta4, kjpindex, index)   
1759          CALL histwrite_p(hist2_id, 'vbeta5', kjit, vbeta5, kjpindex, index)   
1760          CALL histwrite_p(hist2_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
1761          CALL histwrite_p(hist2_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
1762          CALL histwrite_p(hist2_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
1763          IF ( .NOT. hydrol_cwrr ) THEN
1764             CALL histwrite_p(hist2_id, 'rsol', kjit, rsol, kjpindex, index)
1765          ENDIF
1766          IF (.NOT. ( river_routing .AND. nbp_glo .GT. 1) ) THEN
1767          ! if routing not activated output here
1768          ! otherwise output is inside routing
1769               CALL histwrite_p(hist2_id, 'irrigation', kjit, irrigation, kjpindex,index)
1770          ENDIF
1771          CALL histwrite_p(hist2_id, 'snow', kjit, snow, kjpindex, index)
1772          CALL histwrite_p(hist2_id, 'snowage', kjit, snow_age, kjpindex, index)
1773          CALL histwrite_p(hist2_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
1774          CALL histwrite_p(hist2_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
1775          !
1776          IF (  hydrol_cwrr ) THEN
1777             CALL histwrite_p(hist2_id, 'soilindex',  kjit, REAL(njsc, r_std), kjpindex, index)
1778             CALL histwrite_p(hist2_id, 'reinf_slope',  kjit, reinf_slope, kjpindex, index)
1779          ENDIF
1780          !
1781          IF ( ok_co2 ) THEN
1782             CALL histwrite_p(hist2_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)   
1783             CALL histwrite_p(hist2_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
1784             CALL histwrite_p(hist2_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)   
1785          ENDIF
1786          IF ( ok_stomate ) THEN
1787             CALL histwrite_p(hist2_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)   
1788          ENDIF
1789       ELSE
1790          ! Write history file in ALMA format
1791          CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
1792          CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
1793          CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
1794          CALL histwrite_p(hist2_id, 'Qf', kjit, fusion, kjpindex, index)
1795          CALL histwrite_p(hist2_id, 'ESoil', kjit, vevapnu, kjpindex, index)
1796          IF ( river_routing .AND. do_floodplains ) THEN
1797             CALL histwrite_p(hist2_id, 'EWater', kjit, vevapflo, kjpindex, index)
1798             CALL histwrite_p(hist2_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
1799          ENDIF
1800          CALL histwrite_p(hist2_id, 'SWE', kjit, snow, kjpindex, index)
1801          histvar(:)=zero
1802          DO jv=1,nvm
1803             histvar(:) = histvar(:) + transpir(:,jv)
1804          ENDDO
1805          CALL histwrite_p(hist2_id, 'TVeg', kjit, histvar, kjpindex, index)
1806          histvar(:)=zero
1807          DO jv=1,nvm
1808             histvar(:) = histvar(:) + vevapwet(:,jv)
1809          ENDDO
1810          CALL histwrite_p(hist2_id, 'ECanop', kjit, histvar, kjpindex, index)
1811          CALL histwrite_p(hist2_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
1812          CALL histwrite_p(hist2_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
1813          IF ( ok_co2 ) THEN
1814             CALL histwrite_p(hist2_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
1815          ENDIF
1816          IF ( ok_stomate ) THEN
1817             CALL histwrite_p(hist2_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)   
1818          ENDIF
1819       ENDIF ! almaoutput
1820    ENDIF ! hist2_id
1821
1822    !! Change the vegetation fractions if a new map was read in slowproc. This is done
1823    !! after lcchange has been done in stomatelpj
1824        !note by xuhui: according to my reading, putting slowproc_change_frac here is
1825        !logically problematic, because in the end of slowproc, slowproc_veget has cleaned
1826        !carbon content of PFT with veget_max==0. So the new PFT will always
1827        !have 0 Carbon content when started, instead of inheriting C of reduced PFTs
1828    IF (done_stomate_lcchange) THEN
1829       IF (.NOT. use_age_class) THEN
1830           CALL slowproc_change_frac(kjpindex, f_rot_sech, lai, &
1831                                 veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile, fraclut, nwdFraclut)
1832       ENDIF
1833       done_stomate_lcchange=.FALSE.
1834    END IF
1835
1836    !! 14. If it is the last time step, write restart files
1837    IF (ldrestart_write) THEN
1838       CALL sechiba_finalize( &
1839            kjit,     kjpij,  kjpindex, index,   rest_id, &
1840            tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad, &
1841            netco2flux )
1842    END IF
1843
1844  END SUBROUTINE sechiba_main
1845
1846
1847!!  =============================================================================================================================
1848!! SUBROUTINE:    sechiba_finalize
1849!!
1850!>\BRIEF          Finalize all modules by calling their "_finalize" subroutines.
1851!!
1852!! DESCRIPTION:   Finalize all modules by calling their "_finalize" subroutines. These subroutines will write variables to
1853!!                restart file.
1854!!
1855!! \n
1856!_ ==============================================================================================================================
1857
1858  SUBROUTINE sechiba_finalize( &
1859       kjit,     kjpij,  kjpindex, index,   rest_id, &
1860       tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad, &
1861       netco2flux )
1862
1863!! 0.1 Input variables   
1864    INTEGER(i_std), INTENT(in)                               :: kjit              !! Time step number (unitless)
1865    INTEGER(i_std), INTENT(in)                               :: kjpij             !! Total size of the un-compressed grid
1866                                                                                  !! (unitless)
1867    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size - terrestrial pixels only
1868                                                                                  !! (unitless)
1869    INTEGER(i_std),INTENT (in)                               :: rest_id           !! _Restart_ file identifier (unitless)
1870    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index             !! Indices of the pixels on the map.
1871                                                                                  !! Sechiba uses a reduced grid excluding oceans
1872                                                                                  !! ::index contains the indices of the
1873                                                                                  !! terrestrial pixels only! (unitless)
1874    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: tsol_rad           !! Radiative surface temperature
1875                                                                                  !! @tex $(W m^{-2})$ @endtex
1876    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: vevapp             !! Total of evaporation
1877                                                                                  !! @tex $(kg m^{-2} days^{-1})$ @endtex
1878    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: fluxsens           !! Sensible heat flux
1879                                                                                  !! @tex $(W m^{-2})$ @endtex
1880    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: fluxlat            !! Latent heat flux
1881                                                                                  !! @tex $(W m^{-2})$ @endtex
1882    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: tq_cdrag           !! Surface drag coefficient (-)
1883    REAL(r_std),DIMENSION (kjpindex), INTENT (out)          :: netco2flux         !! Sum CO2 flux over PFTs
1884
1885!! 0.2 Local variables
1886    INTEGER(i_std)                                          :: ji, jv             !! Index (unitless)
1887    REAL(r_std), DIMENSION(kjpindex)                        :: histvar            !! Computations for history files (unitless)
1888    CHARACTER(LEN=80)                                       :: var_name           !! To store variables names for I/O (unitless)
1889
1890
1891    !! Write restart file for the next simulation from SECHIBA and other modules
1892
1893    IF (printlev_loc>=3) WRITE (numout,*) 'Start sechiba_finalize for writing restart files'
1894
1895    !! 1. Call diffuco_finalize to write restart files
1896    CALL diffuco_finalize (kjit, kjpindex, rest_id, rstruct, tq_cdrag_pft)
1897   
1898    !! 2. Call energy budget to write restart files
1899    CALL enerbil_finalize (kjit,   kjpindex,    rest_id,            &
1900                           evapot, evapot_corr, temp_sol, temp_sol_pft,  tsol_rad, &
1901                           qsurf,  fluxsens,    fluxlat,  vevapp )
1902   
1903    !! 3. Call hydrology to write restart files
1904    IF ( .NOT. hydrol_cwrr ) THEN
1905       !! 3.1 Call water balance from Choisnel module (2 soil layers) to write restart file
1906         CALL hydrolc_finalize( kjit,      kjpindex,   rest_id,        snow,  &
1907                         snow_age,   snow_nobio, snow_nobio_age, humrel,      &
1908                         vegstress,  qsintveg,   snowrho,        snowtemp,    &
1909                         snowdz,     snowheat,   snowgrain,                   &
1910                         drysoil_frac,           rsol,           shumdiag)
1911       evap_bare_lim(:) = -un
1912       k_litt(:) = huit
1913       shumdiag_perma(:,:)=shumdiag(:,:)
1914    ELSE
1915       !! 3.2 Call water balance from CWRR module (11 soil layers) to write restart file
1916        CALL hydrol_finalize( kjit,           kjpindex, rest_id,  vegstress,  &
1917                        qsintveg,       humrel,         snow,     snow_age, snow_nobio, &
1918                        snow_nobio_age, snowrho,  snowtemp,             &
1919                        snowdz,         snowheat,       &
1920                        snowgrain,  drysoil_frac, evap_bare_lim)
1921!                        snowcap,        snowgrain,  drysoil_frac, evap_bare_lim, evap_bare_lim_pft)
1922       rsol(:) = -un
1923    ENDIF
1924   
1925    !! 4. Call condveg to write surface variables to restart files
1926    CALL condveg_finalize (kjit, kjpindex, rest_id, z0m, z0h, roughheight, roughheight_pft)
1927   
1928    !! 5. Call soil thermodynamic to write restart files
1929    IF (hydrol_cwrr) THEN
1930       CALL thermosoil_finalize (kjit,    kjpindex, rest_id,   gtemp, &
1931            soilcap, soilcap_pft, soilflx, soilflx_pft, lambda_snow, cgrnd_snow, dgrnd_snow)
1932    ELSE
1933       CALL thermosoilc_finalize (kjit,    kjpindex, rest_id,   gtemp, &
1934            soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
1935    END IF
1936
1937    !! 6. Add river routing to restart files 
1938    IF ( river_routing .AND. nbp_glo .GT. 1) THEN
1939       !! 6.1 Call river routing to write restart files
1940       CALL routing_finalize( kjit, kjpindex, rest_id, flood_frac, flood_res )
1941    ELSE
1942       !! 6.2 No routing, set variables to zero
1943       reinfiltration(:) = zero
1944       returnflow(:) = zero
1945!       irrigation(:) = zero
1946       flood_frac(:) = zero
1947       flood_res(:) = zero
1948       IF (do_fullirr) THEN
1949           CALL restput_p (rest_id, 'irrigation', nbp_glo, 1, 1, kjit, irrigation, 'scatter',  nbp_glo, index_g)
1950       ELSE
1951           irrigation(:) = zero
1952       ENDIF
1953    ENDIF
1954   
1955    !! 7. Call slowproc_main to add 'daily' and annual variables to restart file
1956    CALL slowproc_finalize (kjit,       kjpindex,  rest_id,  index,      &
1957                           njsc,       lai,&
1958!!xchen
1959               !!            lai_cl1, lai_cl2, lai_cl3, lai_cl4,&
1960                           height,   veget,      &
1961                           frac_nobio, veget_max, reinf_slope,&
1962                           zz_deep, zz_coef_deep, thawed_humidity, depth_organic_soil, &
1963                           assim_param, frac_age,altmax)
1964
1965    ! Compute global CO2 flux !*
1966    netco2flux(:) = zero
1967    DO jv = 2,nvm
1968      netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*veget_max(:,jv)
1969    ENDDO
1970   
1971    IF (printlev_loc>=3) WRITE (numout,*) 'sechiba_finalize done'
1972   
1973  END SUBROUTINE sechiba_finalize
1974
1975 
1976!! ==============================================================================================================================\n
1977!! SUBROUTINE   : sechiba_init
1978!!
1979!>\BRIEF        Dynamic allocation of the variables, the dimensions of the
1980!! variables are determined by user-specified settings.
1981!!
1982!! DESCRIPTION  : The domain size (:: kjpindex) is used to allocate the correct
1983!! dimensions to all variables in sechiba. Depending on the variable, its
1984!! dimensions are also determined by the number of PFT's (::nvm), number of
1985!! soil types (::nstm), number of non-vegetative surface types (::nnobio),
1986!! number of soil levels (::ngrnd), number of soil layers in the hydrological
1987!! model (i.e. cwrr) (::nslm). Values for these variables are set in
1988!! constantes_soil.f90 and constantes_veg.f90.\n
1989!!
1990!! Memory is allocated for all Sechiba variables and new indexing tables
1991!! are build making use of both (::kjpij) and (::kjpindex). New indexing tables
1992!! are needed because a single pixel can contain several PFTs, soil types, etc.
1993!! The new indexing tables have separate indices for the different
1994!! PFTs, soil types, etc.\n
1995!!
1996!! RECENT CHANGE(S): None
1997!!
1998!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
1999!! variables. However, the routine allocates memory and builds new indexing
2000!! variables for later use.
2001!!
2002!! REFERENCE(S) : None
2003!!
2004!! FLOWCHART    : None
2005!! \n
2006!_ ================================================================================================================================
2007
2008  SUBROUTINE sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)
2009
2010!! 0.1 Input variables
2011 
2012    INTEGER(i_std), INTENT (in)                         :: kjit               !! Time step number (unitless)
2013    INTEGER(i_std), INTENT (in)                         :: kjpij              !! Total size of the un-compressed grid (unitless)
2014    INTEGER(i_std), INTENT (in)                         :: kjpindex           !! Domain size - terrestrial pixels only (unitless)
2015    INTEGER(i_std), INTENT (in)                         :: rest_id            !! _Restart_ file identifier (unitless)
2016    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: index              !! Indeces of the points on the map (unitless)
2017    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)     :: lalo               !! Geographical coordinates (latitude,longitude)
2018                                                                              !! for pixels (degrees)
2019!! 0.2 Output variables
2020
2021!! 0.3 Modified variables
2022
2023!! 0.4 Local variables
2024
2025    INTEGER(i_std)                                      :: ier                !! Check errors in memory allocation (unitless)
2026    INTEGER(i_std)                                      :: ji, jv             !! Indeces (unitless)
2027    REAL(r_std), DIMENSION (ngrnd)                      :: dz_tmp             !! Dummy variable
2028!_ ==============================================================================================================================
2029
2030!! 1. Initialize variables
2031   
2032    ! Dynamic allocation with user-specified dimensions on first call
2033    IF (l_first_sechiba) THEN
2034       l_first_sechiba=.FALSE.
2035    ELSE
2036       CALL ipslerr_p(3,'sechiba_init',' l_first_sechiba false . we stop ','','')
2037    ENDIF
2038
2039    !! Initialize local printlev
2040    printlev_loc=get_printlev('sechiba')
2041   
2042
2043    !! 1.1 Initialize 3D vegetation indexation table
2044    ALLOCATE (indexveg(kjpindex*nvm),stat=ier)
2045    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexveg','','')
2046
2047    ALLOCATE (indexlai(kjpindex*(nlai+1)),stat=ier)
2048    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlai','','')
2049
2050    ALLOCATE (indexlai0(kjpindex*(nlai)),stat=ier)
2051    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlai0','','')
2052
2053    ALLOCATE (indexsoil(kjpindex*nstm),stat=ier)
2054    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsoil','','')
2055
2056    ALLOCATE (indexnobio(kjpindex*nnobio),stat=ier)
2057    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnobio','','')
2058
2059    ALLOCATE (indexgrnd(kjpindex*ngrnd),stat=ier)
2060    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexgrnd','','')
2061
2062    ALLOCATE (indexsnow(kjpindex*nsnow),stat=ier)
2063    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsnow','','')
2064
2065    ALLOCATE (indexlayer(kjpindex*nslm),stat=ier)
2066    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlayer','','')
2067
2068    ALLOCATE (indexnslm(kjpindex*nslm),stat=ier)
2069    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnslm','','')
2070
2071    ALLOCATE (indexalb(kjpindex*2),stat=ier)
2072    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexalb','','')
2073
2074    !! 1.2  Initialize 1D array allocation with restartable value
2075    ALLOCATE (flood_res(kjpindex),stat=ier)
2076    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for flood_res','','')
2077    flood_res(:) = undef_sechiba
2078
2079    ALLOCATE (flood_frac(kjpindex),stat=ier)
2080    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for kjpindex','','')
2081    flood_frac(:) = undef_sechiba
2082
2083    ALLOCATE (snow(kjpindex),stat=ier)
2084    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow','','')
2085    snow(:) = undef_sechiba
2086
2087    ALLOCATE (snow_age(kjpindex),stat=ier)
2088    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_age','','')
2089    snow_age(:) = undef_sechiba
2090
2091    ALLOCATE (drysoil_frac(kjpindex),stat=ier)
2092    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drysoil_frac','','')
2093
2094    ALLOCATE (rsol(kjpindex),stat=ier)
2095    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rsol','','')
2096
2097    ALLOCATE (evap_bare_lim(kjpindex),stat=ier)
2098    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim','','')
2099
2100!    ALLOCATE (evap_bare_lim_pft(kjpindex,nvm),stat=ier)
2101!    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim_pft','','')
2102
2103    ALLOCATE (soil_deficit(kjpindex,nvm),stat=ier)
2104    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soil_deficit','','')
2105
2106    ALLOCATE (evapot(kjpindex),stat=ier)
2107    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot','','')
2108    evapot(:) = undef_sechiba
2109
2110    ALLOCATE (evapot_corr(kjpindex),stat=ier)
2111    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot_corr','','')
2112
2113    ALLOCATE (humrel(kjpindex,nvm),stat=ier)
2114    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for humrel','','')
2115    humrel(:,:) = undef_sechiba
2116
2117    ALLOCATE (vegstress(kjpindex,nvm),stat=ier)
2118    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vegstress','','')
2119    vegstress(:,:) = undef_sechiba
2120
2121    ALLOCATE (vegstress_old(kjpindex,nvm),stat=ier)
2122    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vegstress_old','','')
2123    vegstress_old(:,:) = 1
2124
2125    ALLOCATE (njsc(kjpindex),stat=ier)
2126    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for njsc','','')
2127    njsc(:)= undef_int
2128
2129    ALLOCATE (soiltile(kjpindex,nstm),stat=ier)
2130    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soiltile','','')
2131
2132    ALLOCATE (is_crop_soil(nstm),stat=ier)
2133    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for is_crop_soil','','')
2134
2135    ALLOCATE (fraclut(kjpindex,nlut),stat=ier)
2136    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fraclut','','')
2137
2138    ALLOCATE (nwdFraclut(kjpindex,nlut),stat=ier)
2139    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for nwdFraclut','','')
2140
2141    ALLOCATE (reinf_slope(kjpindex),stat=ier)
2142    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope','','')
2143
2144    ALLOCATE (vbeta1(kjpindex),stat=ier)
2145    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta1','','')
2146
2147    ALLOCATE (vbeta4(kjpindex),stat=ier)
2148    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta4','','')
2149
2150    ALLOCATE (vbeta4_pft(kjpindex,nvm),stat=ier)
2151    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta4_pft','','')
2152
2153    ALLOCATE (vbeta5(kjpindex),stat=ier)
2154    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta5','','')
2155
2156    ALLOCATE (soilcap(kjpindex),stat=ier)
2157    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilcap','','')
2158
2159    ALLOCATE (soilcap_pft(kjpindex,nvm),stat=ier)
2160    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilcap_pft','','')
2161
2162    ALLOCATE (soilflx(kjpindex),stat=ier)
2163    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilflx','','')
2164
2165    ALLOCATE (soilflx_pft(kjpindex,nvm),stat=ier)
2166    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilflx_pft','','')
2167
2168    ALLOCATE (temp_sol(kjpindex),stat=ier)
2169    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol','','')
2170    temp_sol(:) = undef_sechiba
2171
2172    ALLOCATE (temp_sol_pft(kjpindex, nvm),stat=ier)
2173    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol_pft','','')
2174    temp_sol_pft(:,:) = undef_sechiba
2175
2176    ALLOCATE (temp_sol_new_pft(kjpindex, nvm),stat=ier)
2177    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol_new_pft','','')
2178    temp_sol_new_pft(:,:) = undef_sechiba
2179
2180
2181    ALLOCATE (qsurf(kjpindex),stat=ier)
2182    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsurf','','')
2183    qsurf(:) = undef_sechiba
2184
2185    !! 1.3 Initialize 2D array allocation with restartable value
2186    ALLOCATE (qsintveg(kjpindex,nvm),stat=ier)
2187    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintveg','','')
2188    qsintveg(:,:) = undef_sechiba
2189
2190    ALLOCATE (tq_cdrag_pft(kjpindex,nvm),stat=ier)
2191    IF (ier.NE.0) THEN
2192       WRITE (numout,*) ' error in tq_cdrag_pft allocation. We stop. We need kjpindex x nvm words = ',&
2193            & kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
2194       STOP 'sechiba_init'
2195    END IF
2196    ALLOCATE (vbeta_pft(kjpindex,nvm),stat=ier)
2197    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta_pft','','')
2198
2199    ALLOCATE (vbeta2(kjpindex,nvm),stat=ier)
2200    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta2','','')
2201
2202    ALLOCATE (vbeta3(kjpindex,nvm),stat=ier)
2203    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3','','')
2204
2205    ALLOCATE (vbeta3pot(kjpindex,nvm),stat=ier)
2206    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3pot','','')
2207
2208    ALLOCATE (gsmean(kjpindex,nvm),stat=ier)
2209    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gsmean','','')
2210
2211    ALLOCATE (cimean(kjpindex,nvm),stat=ier)
2212    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cimean','','')
2213
2214    ALLOCATE (gpp(kjpindex,nvm),stat=ier)
2215    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gpp','','')
2216    gpp(:,:) = undef_sechiba
2217
2218 
2219    ALLOCATE (temp_growth(kjpindex),stat=ier) 
2220    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_growth','','')
2221    temp_growth(:) = undef_sechiba 
2222
2223    ALLOCATE (veget(kjpindex,nvm),stat=ier)
2224    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget','','')
2225    veget(:,:)=undef_sechiba
2226
2227    ALLOCATE (veget_max(kjpindex,nvm),stat=ier)
2228    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget_max','','')
2229
2230    ALLOCATE (tot_bare_soil(kjpindex),stat=ier)
2231    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_bare_soil','','')
2232
2233    ALLOCATE (lai(kjpindex,nvm),stat=ier)
2234    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lai','','')
2235    lai(:,:)=undef_sechiba
2236
2237!!@chy: lai_cl1/2/3/4 has to be allocated here. @xchen
2238    ALLOCATE (lai_cl1(kjpindex,nvm),stat=ier)                                       
2239    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lai_cl1','','') 
2240    lai_cl1(:,:)=undef_sechiba
2241    ALLOCATE (lai_cl2(kjpindex,nvm),stat=ier)                                   
2242    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lai_cl2','','') 
2243    lai_cl2(:,:)=undef_sechiba
2244    ALLOCATE (lai_cl3(kjpindex,nvm),stat=ier)                                   
2245    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lai_cl3','','') 
2246    lai_cl3(:,:)=undef_sechiba 
2247    ALLOCATE (lai_cl4(kjpindex,nvm),stat=ier)                                   
2248    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lai_cl4','','') 
2249    lai_cl4(:,:)=undef_sechiba 
2250    ALLOCATE (VPD(kjpindex),stat=ier)                                   
2251    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for VPD','','')
2252    VPD(:)=undef_sechiba 
2253    ALLOCATE (ave_VPD(kjpindex),stat=ier)                                           
2254    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ave_VPD','','') 
2255    ave_VPD(:)=undef_sechiba   
2256    ALLOCATE (laitop(kjpindex,nvm),stat=ier)                                   
2257    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for laitop','','')
2258    laitop(:,:)=undef_sechiba                                                   
2259    ALLOCATE (laidow(kjpindex,nvm),stat=ier)                                     
2260    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for laidow','','')
2261    laidow(:,:)=undef_sechiba 
2262    ALLOCATE (vcmax_cl1(kjpindex,nvm),stat=ier)                                   
2263    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vcmax_cl1','','')
2264    vcmax_cl1(:,:)=undef_sechiba 
2265    ALLOCATE (vcmax_cl2(kjpindex,nvm),stat=ier)                                 
2266    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vcmax_cl2','','')
2267    vcmax_cl2(:,:)=undef_sechiba 
2268    ALLOCATE (vcmax_cl3(kjpindex,nvm),stat=ier)                                 
2269    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vcmax_cl3','','')
2270    vcmax_cl3(:,:)=undef_sechiba 
2271    ALLOCATE (vcmax_cl4(kjpindex,nvm),stat=ier)                                 
2272    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vcmax_cl4','','')
2273    vcmax_cl4(:,:)=undef_sechiba 
2274
2275    ALLOCATE (frac_age(kjpindex,nvm,nleafages),stat=ier)
2276    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_age','','')
2277    frac_age(:,:,:)=undef_sechiba
2278
2279    ALLOCATE (height(kjpindex,nvm),stat=ier)
2280    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for height','','')
2281    height(:,:)=undef_sechiba
2282
2283    ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
2284    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_nobio','','')
2285    frac_nobio(:,:) = undef_sechiba
2286
2287    ALLOCATE (snow_nobio(kjpindex,nnobio),stat=ier)
2288    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio','','')
2289    snow_nobio(:,:) = undef_sechiba
2290
2291    ALLOCATE (snow_nobio_age(kjpindex,nnobio),stat=ier)
2292    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio_age','','')
2293    snow_nobio_age(:,:) = undef_sechiba
2294
2295    ALLOCATE (assim_param(kjpindex,nvm,npco2),stat=ier)
2296    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for assim_param','','')
2297
2298    !! 1.4 Initialize 1D array allocation
2299! xuhui: +
2300    ALLOCATE(f_rot_sech(kjpindex),stat=ier)
2301    IF (ier.NE.0) THEN
2302       WRITE (numout,*) ' error in f_rot_sech allocation. We stop. we need kjpindex words = ',kjpindex
2303       STOP 'sechiba_init'
2304    END IF
2305    f_rot_sech(:) = .FALSE.
2306   
2307! xuhui: -
2308!pss:+
2309!    ALLOCATE (fwet_out(kjpindex),stat=ier)
2310!    IF (ier.NE.0) THEN
2311!       WRITE (numout,*) ' error in fwet_out allocation. We stop. we need kjpindex words = ',kjpindex
2312!       STOP 'sechiba_init'
2313!    END IF
2314!    fwet_out(:) = undef_sechiba
2315!pss:-
2316!pss:+
2317!    ALLOCATE (drunoff_tot(kjpindex),stat=ier)
2318!    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drunoff_tot','','')
2319!    drunoff_tot(:) = zero
2320!pss:-
2321
2322    ALLOCATE (vevapflo(kjpindex),stat=ier)
2323    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapflo','','')
2324    vevapflo(:)=zero
2325
2326    ALLOCATE (vevapsno(kjpindex),stat=ier)
2327    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapsno','','')
2328
2329    ALLOCATE (vevapnu(kjpindex),stat=ier)
2330    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapnu','','')
2331
2332    ALLOCATE (vevapnu_pft(kjpindex,nvm),stat=ier)
2333    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapnu_pft','','')
2334
2335    ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
2336    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for totfrac_nobio','','')
2337
2338    ALLOCATE (floodout(kjpindex),stat=ier)
2339    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for floodout','','')
2340
2341    ALLOCATE (runoff(kjpindex),stat=ier)
2342    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for runoff','','')
2343
2344    ALLOCATE (drainage(kjpindex),stat=ier)
2345    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drainage','','')
2346
2347    ALLOCATE (returnflow(kjpindex),stat=ier)
2348    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for returnflow','','')
2349    returnflow(:) = zero
2350
2351    ALLOCATE (reinfiltration(kjpindex),stat=ier)
2352    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinfiltration','','')
2353    reinfiltration(:) = zero
2354
2355    ALLOCATE (irrigation(kjpindex),stat=ier)
2356    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for irrigation','','')
2357    irrigation(:) = zero
2358
2359    IF ( do_fullirr ) THEN
2360        ALLOCATE (irrig_frac(kjpindex),stat=ier)
2361        IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for irrig_frac','','')
2362        irrig_frac(:) = 1 ! percentage of cropland irrigated
2363          !! put here a irrigation map input routine, xuhui
2364          !!
2365
2366!              ALLOCATE (tot_vegfrac_crop(kjpindex),stat=ier)
2367!              IF (ier.NE.0) THEN
2368!                 WRITE (numout,*) ' error in tot_vegfrac_crop allocation. We stop. We need kjpindex words = ',kjpindex
2369!                 STOP 'sechiba_main init'
2370!              END IF
2371!              DO jv  = 2, nvm
2372!                 IF ( (jv /= ibare_sechiba) .AND. .NOT.(natural(jv)) ) THEN
2373!                    tot_vegfrac_crop(:) = tot_vegfrac_crop(:) + veget_max(:,jv)
2374!                 END IF
2375!              END DO
2376    ENDIF
2377
2378    ALLOCATE (z0h(kjpindex),stat=ier)
2379    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0h','','')
2380
2381    ALLOCATE (z0m(kjpindex),stat=ier)
2382    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0m','','')
2383
2384    ALLOCATE (roughheight(kjpindex),stat=ier)
2385    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for roughheight','','')
2386
2387    ALLOCATE (roughheight_pft(kjpindex,nvm),stat=ier) 
2388    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for roughheight_pft','','')
2389
2390    ALLOCATE (emis(kjpindex),stat=ier)
2391    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for emis','','')
2392
2393    ALLOCATE (tot_melt(kjpindex),stat=ier)
2394    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_melt','','')
2395
2396    ALLOCATE (vbeta(kjpindex),stat=ier)
2397    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta','','')
2398
2399    ALLOCATE (fusion(kjpindex),stat=ier)
2400    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fusion','','')
2401
2402    ALLOCATE (rau(kjpindex),stat=ier)
2403    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rau','','')
2404
2405    ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
2406    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for deadleaf_cover','','')
2407
2408    ALLOCATE (stempdiag(kjpindex, nslm),stat=ier)
2409    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for stempdiag','','')
2410
2411    ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
2412    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for co2_flux','','')
2413    co2_flux(:,:)=zero
2414
2415    ALLOCATE (shumdiag(kjpindex,nslm),stat=ier)
2416    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag','','')
2417   
2418    ALLOCATE (shumdiag_perma(kjpindex,nslm),stat=ier)
2419    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag_perma','','')
2420
2421    ALLOCATE (litterhumdiag(kjpindex),stat=ier)
2422    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for litterhumdiag','','')
2423!gmjc top 5 layer grassland soil moisture for grazing
2424    ALLOCATE (tmc_topgrass(kjpindex),stat=ier)
2425    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tmc_topgrass','','')
2426!end gmjc
2427
2428    ALLOCATE (humcste_use(kjpindex,nvm),stat=ier)
2429    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for humcste_use','','')
2430
2431    ALLOCATE (altmax(kjpindex,nvm),stat=ier)
2432    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for altmax','','')
2433
2434    ALLOCATE (ptnlev1(kjpindex),stat=ier)
2435    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ptnlev1','','')
2436
2437    ALLOCATE (k_litt(kjpindex),stat=ier)
2438    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for k_litt','','')
2439
2440    !! 1.5 Initialize 2D array allocation
2441    ALLOCATE (vevapwet(kjpindex,nvm),stat=ier)
2442    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapwet','','')
2443    vevapwet(:,:)=undef_sechiba
2444
2445    ALLOCATE (transpir(kjpindex,nvm),stat=ier)
2446    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpir','','')
2447
2448    ALLOCATE (transpot(kjpindex,nvm),stat=ier)
2449    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpot','','')
2450
2451    ALLOCATE (qsintmax(kjpindex,nvm),stat=ier)
2452    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintmax','','')
2453
2454    ALLOCATE (rveget(kjpindex,nvm),stat=ier)
2455    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rveget','','')
2456
2457    ALLOCATE (rstruct(kjpindex,nvm),stat=ier)
2458    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rstruct','','')
2459
2460    ALLOCATE (pgflux(kjpindex),stat=ier)
2461    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for pgflux','','')
2462    pgflux(:)= 0.0
2463
2464    ALLOCATE (soilflxresid(kjpindex),stat=ier)
2465    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilflxresid','','')
2466    soilflxresid = zero
2467
2468    ALLOCATE (cgrnd_snow(kjpindex,nsnow),stat=ier)
2469    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cgrnd_snow','','')
2470    cgrnd_snow(:,:) = 0
2471
2472    ALLOCATE (dgrnd_snow(kjpindex,nsnow),stat=ier)
2473    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for dgrnd_snow','','')
2474    dgrnd_snow(:,:) = 0
2475
2476    ALLOCATE (lambda_snow(kjpindex),stat=ier)
2477    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lambda_snow','','')
2478    lambda_snow(:) = 0
2479
2480    ALLOCATE (frac_snow_veg(kjpindex),stat=ier)
2481    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_snow_veg','','')
2482
2483    ALLOCATE (frac_snow_nobio(kjpindex,nnobio),stat=ier)
2484    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_snow_nobio','','')
2485
2486    ALLOCATE (snowrho(kjpindex,nsnow),stat=ier)
2487    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowrho','','')
2488    snowrho = zero
2489
2490    ALLOCATE (snowheat(kjpindex,nsnow),stat=ier)
2491    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowheat','','')
2492
2493    ALLOCATE (snowgrain(kjpindex,nsnow),stat=ier)
2494    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowgrain','','')
2495
2496    ALLOCATE (snowtemp(kjpindex,nsnow),stat=ier)
2497    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowtemp','','')
2498
2499    ALLOCATE (snowdz(kjpindex,nsnow),stat=ier)
2500    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowdz','','')
2501    snowdz = zero
2502
2503    ALLOCATE (temp_sol_add(kjpindex),stat=ier)
2504    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol_add','','')
2505    temp_sol_add(:) = 0.0
2506
2507    ALLOCATE (gtemp(kjpindex),stat=ier)
2508    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gtemp','','')
2509
2510    !allocate arrays needed for permafrost calculations
2511    ALLOCATE(tdeep(kjpindex,ndeep,nvm),stat=ier)
2512    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tdeep','','')
2513    tdeep(:,:,:) = 250.
2514
2515    ALLOCATE(hsdeep(kjpindex,ndeep,nvm),stat=ier)
2516    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for hsdeep','','')
2517    hsdeep(:,:,:) = 1.0
2518
2519    ALLOCATE(heat_Zimov(kjpindex,ndeep,nvm),stat=ier)
2520    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for heat_Zimov','','')
2521    heat_Zimov(:,:,:) = zero
2522
2523  ! 1d arrays (xy)
2524    ALLOCATE(sfluxCH4_deep(kjpindex),stat=ier)
2525    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for sfluxCH4_deep','','')
2526    ALLOCATE(sfluxCO2_deep(kjpindex),stat=ier)
2527    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for sfluxCO2_deep','','')
2528    ALLOCATE(thawed_humidity(kjpindex),stat=ier)
2529    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for thawed_humidity','','')
2530    ALLOCATE(depth_organic_soil(kjpindex),stat=ier)
2531    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for depth_organic_soil','','')
2532    ! 1d arrays (ndeep)
2533    ALLOCATE(zz_deep(ndeep),stat=ier)
2534    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for zz_deep','','')
2535    ALLOCATE(zz_coef_deep(ndeep),stat=ier)
2536    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for zz_coef_deep','','')
2537    ALLOCATE(soilc_total(kjpindex,ndeep,nvm),stat=ier)
2538    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilc_total','','')
2539
2540    ALLOCATE (mc_layh(kjpindex, nslm),stat=ier)
2541    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mc_layh','','')
2542
2543    ALLOCATE (mcl_layh(kjpindex, nslm),stat=ier)
2544    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcl_layh','','')
2545
2546    ALLOCATE (soilmoist(kjpindex, nslm),stat=ier)
2547    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilmoist','','')
2548
2549    ALLOCATE (mc_layh_s(kjpindex, nslm, nstm),stat=ier)
2550    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mc_layh_s','','')
2551
2552    ALLOCATE (mcl_layh_s(kjpindex, nslm, nstm),stat=ier)
2553    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcl_layh_s','','')
2554
2555!    ALLOCATE (mc_layh_pft(kjpindex, nslm, nvm),stat=ier)
2556!    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mc_layh_pft','','')
2557
2558!    ALLOCATE (mcl_layh_pft(kjpindex, nslm, nvm),stat=ier)
2559!    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcl_layh_pft','','')
2560
2561!    ALLOCATE (soilmoist_pft(kjpindex, nslm, nvm),stat=ier)
2562!    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilmoist_pft','','')
2563
2564    !! 1.6 Initialize indexing table for the vegetation fields.
2565    ! In SECHIBA we work on reduced grids but to store in the full 3D filed vegetation variable
2566    ! we need another index table : indexveg, indexsoil, indexnobio and indexgrnd
2567    DO ji = 1, kjpindex
2568       !
2569       DO jv = 1, nlai
2570          indexlai0((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2571       ENDDO
2572       !
2573       DO jv = 1, nlai+1
2574          indexlai((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2575       ENDDO
2576       !
2577       DO jv = 1, nvm
2578          indexveg((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2579       ENDDO
2580       !     
2581       DO jv = 1, nstm
2582          indexsoil((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2583       ENDDO
2584       !     
2585       DO jv = 1, nnobio
2586          indexnobio((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2587       ENDDO
2588       !
2589       DO jv = 1, ngrnd
2590          indexgrnd((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2591       ENDDO
2592       !
2593       DO jv = 1, nsnow
2594          indexsnow((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2595       ENDDO
2596
2597       DO jv = 1, nslm
2598          indexnslm((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2599       ENDDO
2600
2601       DO jv = 1, nslm
2602          indexlayer((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2603       ENDDO
2604       !
2605       DO jv = 1, 2
2606          indexalb((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
2607       ENDDO
2608       !
2609    ENDDO
2610    !
2611    zz_deep      = znt  ! Vertical depth of the nodes for Thermodynamics
2612    zz_coef_deep = zlt  ! Vertical depth of the Layers for Thermodynamics
2613    dz_tmp       = dlt  ! Delta of the Layers for Thermodynamics
2614    !
2615
2616!! 2. Read the default value that will be put into variable which are not in the restart file
2617    CALL ioget_expval(val_exp)
2618   
2619    IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_init done '
2620
2621  END SUBROUTINE sechiba_init
2622 
2623
2624!! ==============================================================================================================================\n
2625!! SUBROUTINE   : sechiba_clear
2626!!
2627!>\BRIEF        Deallocate memory of sechiba's variables
2628!!
2629!! DESCRIPTION  : None
2630!!
2631!! RECENT CHANGE(S): None
2632!!
2633!! MAIN OUTPUT VARIABLE(S): None
2634!!
2635!! REFERENCE(S) : None
2636!!
2637!! FLOWCHART    : None
2638!! \n
2639!_ ================================================================================================================================
2640
2641  SUBROUTINE sechiba_clear()
2642
2643!! 1. Initialize first run
2644
2645    l_first_sechiba=.TRUE.
2646
2647!! 2. Deallocate dynamic variables of sechiba
2648
2649    IF ( ALLOCATED (indexveg)) DEALLOCATE (indexveg)
2650    IF ( ALLOCATED (indexlai0)) DEALLOCATE (indexlai0)
2651    IF ( ALLOCATED (indexlai)) DEALLOCATE (indexlai)
2652    IF ( ALLOCATED (indexsoil)) DEALLOCATE (indexsoil)
2653    IF ( ALLOCATED (indexnobio)) DEALLOCATE (indexnobio)
2654    IF ( ALLOCATED (indexsnow)) DEALLOCATE (indexsnow)
2655    IF ( ALLOCATED (indexgrnd)) DEALLOCATE (indexgrnd)
2656    IF ( ALLOCATED (indexlayer)) DEALLOCATE (indexlayer)
2657    IF ( ALLOCATED (indexnslm)) DEALLOCATE (indexnslm)
2658    IF ( ALLOCATED (indexalb)) DEALLOCATE (indexalb)
2659    IF ( ALLOCATED (flood_res)) DEALLOCATE (flood_res)
2660    IF ( ALLOCATED (flood_frac)) DEALLOCATE (flood_frac)
2661    IF ( ALLOCATED (snow)) DEALLOCATE (snow)
2662    IF ( ALLOCATED (snow_age)) DEALLOCATE (snow_age)
2663    IF ( ALLOCATED (drysoil_frac)) DEALLOCATE (drysoil_frac)
2664    IF ( ALLOCATED (rsol)) DEALLOCATE (rsol)
2665    IF ( ALLOCATED (evap_bare_lim)) DEALLOCATE (evap_bare_lim)
2666!    IF ( ALLOCATED (evap_bare_lim_pft)) DEALLOCATE (evap_bare_lim_pft)
2667    IF ( ALLOCATED (soil_deficit)) DEALLOCATE (soil_deficit)
2668    IF ( ALLOCATED (evapot)) DEALLOCATE (evapot)
2669    IF ( ALLOCATED (evapot_corr)) DEALLOCATE (evapot_corr)
2670    IF ( ALLOCATED (humrel)) DEALLOCATE (humrel)
2671    IF ( ALLOCATED (vegstress)) DEALLOCATE (vegstress)
2672    IF ( ALLOCATED (vegstress_old)) DEALLOCATE (vegstress_old)
2673    IF ( ALLOCATED (soiltile)) DEALLOCATE (soiltile)
2674    IF ( ALLOCATED (is_crop_soil)) DEALLOCATE (is_crop_soil)
2675    IF ( ALLOCATED (fraclut)) DEALLOCATE (fraclut)
2676    IF ( ALLOCATED (nwdFraclut)) DEALLOCATE (nwdFraclut)
2677    IF ( ALLOCATED (njsc)) DEALLOCATE (njsc)
2678    IF ( ALLOCATED (reinf_slope)) DEALLOCATE (reinf_slope)
2679    IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1)
2680    IF ( ALLOCATED (vbeta_pft)) DEALLOCATE (vbeta_pft)
2681    IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4)
2682    IF ( ALLOCATED (vbeta4_pft)) DEALLOCATE (vbeta4_pft)
2683    IF ( ALLOCATED (vbeta5)) DEALLOCATE (vbeta5)
2684    IF ( ALLOCATED (soilcap)) DEALLOCATE (soilcap)
2685    IF ( ALLOCATED (soilcap_pft)) DEALLOCATE (soilcap_pft)
2686    IF ( ALLOCATED (soilflx)) DEALLOCATE (soilflx)
2687    IF ( ALLOCATED (soilflx_pft)) DEALLOCATE (soilflx_pft)
2688!    IF ( ALLOCATED (snowcap)) DEALLOCATE (snowcap)
2689!    IF ( ALLOCATED (snowflx)) DEALLOCATE (snowflx)
2690    IF ( ALLOCATED (temp_sol)) DEALLOCATE (temp_sol)
2691    IF ( ALLOCATED (temp_sol_pft)) DEALLOCATE (temp_sol_pft)
2692    IF ( ALLOCATED (qsurf)) DEALLOCATE (qsurf)
2693    IF ( ALLOCATED (qsintveg)) DEALLOCATE (qsintveg)
2694    IF ( ALLOCATED (vbeta2))  DEALLOCATE (vbeta2)
2695    IF ( ALLOCATED (vbeta3)) DEALLOCATE (vbeta3)
2696    IF ( ALLOCATED (vbeta3pot)) DEALLOCATE (vbeta3pot)
2697    IF ( ALLOCATED (tq_cdrag_pft))  DEALLOCATE (tq_cdrag_pft)
2698    IF ( ALLOCATED (gsmean)) DEALLOCATE (gsmean)
2699    IF ( ALLOCATED (cimean)) DEALLOCATE (cimean)
2700    IF ( ALLOCATED (gpp)) DEALLOCATE (gpp)
2701    IF ( ALLOCATED (temp_growth)) DEALLOCATE (temp_growth) 
2702    IF ( ALLOCATED (veget)) DEALLOCATE (veget)
2703    IF ( ALLOCATED (veget_max)) DEALLOCATE (veget_max)
2704    IF ( ALLOCATED (tot_bare_soil)) DEALLOCATE (tot_bare_soil)
2705    IF ( ALLOCATED (lai)) DEALLOCATE (lai)
2706!@chy: these lines below are needed.
2707!!@xchen
2708    IF ( ALLOCATED (lai_cl1)) DEALLOCATE (lai_cl1)                                     
2709    IF ( ALLOCATED (lai_cl2)) DEALLOCATE (lai_cl2)                             
2710    IF ( ALLOCATED (lai_cl3)) DEALLOCATE (lai_cl3)                             
2711    IF ( ALLOCATED (lai_cl4)) DEALLOCATE (lai_cl4)   
2712    IF ( ALLOCATED (VPD)) DEALLOCATE (VPD)   
2713    IF ( ALLOCATED (ave_VPD)) DEALLOCATE (ave_VPD)                                                         
2714    IF ( ALLOCATED (laitop)) DEALLOCATE (laitop)                               
2715    IF ( ALLOCATED (laidow)) DEALLOCATE (laidow) 
2716    IF ( ALLOCATED (vcmax_cl1)) DEALLOCATE (vcmax_cl1)                             
2717    IF ( ALLOCATED (vcmax_cl2)) DEALLOCATE (vcmax_cl2)                         
2718    IF ( ALLOCATED (vcmax_cl3)) DEALLOCATE (vcmax_cl3)                         
2719    IF ( ALLOCATED (vcmax_cl4)) DEALLOCATE (vcmax_cl4)                         
2720
2721
2722    IF ( ALLOCATED (frac_age)) DEALLOCATE (frac_age)
2723    IF ( ALLOCATED (height)) DEALLOCATE (height)
2724    IF ( ALLOCATED (roughheight)) DEALLOCATE (roughheight)
2725    IF ( ALLOCATED (roughheight_pft)) DEALLOCATE (roughheight_pft)
2726    IF ( ALLOCATED (frac_nobio)) DEALLOCATE (frac_nobio)
2727    IF ( ALLOCATED (snow_nobio)) DEALLOCATE (snow_nobio)
2728    IF ( ALLOCATED (snow_nobio_age)) DEALLOCATE (snow_nobio_age)
2729    IF ( ALLOCATED (assim_param)) DEALLOCATE (assim_param)
2730    IF ( ALLOCATED (vevapflo)) DEALLOCATE (vevapflo)
2731    IF ( ALLOCATED (vevapsno)) DEALLOCATE (vevapsno)
2732    IF ( ALLOCATED (vevapnu)) DEALLOCATE (vevapnu)
2733    IF ( ALLOCATED (vevapnu_pft)) DEALLOCATE (vevapnu_pft)
2734    IF ( ALLOCATED (totfrac_nobio)) DEALLOCATE (totfrac_nobio)
2735    IF ( ALLOCATED (floodout)) DEALLOCATE (floodout)
2736    IF ( ALLOCATED (runoff)) DEALLOCATE (runoff)
2737    IF ( ALLOCATED (drainage)) DEALLOCATE (drainage)
2738    IF ( ALLOCATED (reinfiltration)) DEALLOCATE (reinfiltration)
2739    IF ( ALLOCATED (irrigation)) DEALLOCATE (irrigation)
2740    IF ( ALLOCATED (irrig_frac)) DEALLOCATE (irrig_frac)
2741    IF ( ALLOCATED (tot_vegfrac_crop)) DEALLOCATE (tot_vegfrac_crop)
2742    IF ( ALLOCATED (tot_melt)) DEALLOCATE (tot_melt)
2743    IF ( ALLOCATED (vbeta)) DEALLOCATE (vbeta)
2744    IF ( ALLOCATED (fusion)) DEALLOCATE (fusion)
2745    IF ( ALLOCATED (rau)) DEALLOCATE (rau)
2746    IF ( ALLOCATED (deadleaf_cover)) DEALLOCATE (deadleaf_cover)
2747    IF ( ALLOCATED (stempdiag)) DEALLOCATE (stempdiag)
2748    IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
2749    IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
2750    IF ( ALLOCATED (shumdiag_perma)) DEALLOCATE (shumdiag_perma)
2751    IF ( ALLOCATED (litterhumdiag)) DEALLOCATE (litterhumdiag)
2752!gmjc top 5 layer grassland soil moisture for grazing
2753    IF ( ALLOCATED (tmc_topgrass)) DEALLOCATE (tmc_topgrass)
2754!end gmjc
2755    IF ( ALLOCATED (humcste_use)) DEALLOCATE (humcste_use)
2756    IF ( ALLOCATED (altmax)) DEALLOCATE (altmax)
2757    IF ( ALLOCATED (ptnlev1)) DEALLOCATE (ptnlev1)
2758    IF ( ALLOCATED (k_litt)) DEALLOCATE (k_litt)
2759    IF ( ALLOCATED (vevapwet)) DEALLOCATE (vevapwet)
2760    IF ( ALLOCATED (transpir)) DEALLOCATE (transpir)
2761    IF ( ALLOCATED (transpot)) DEALLOCATE (transpot)
2762    IF ( ALLOCATED (qsintmax)) DEALLOCATE (qsintmax)
2763    IF ( ALLOCATED (rveget)) DEALLOCATE (rveget)
2764    IF ( ALLOCATED (rstruct)) DEALLOCATE (rstruct)
2765    IF ( ALLOCATED (frac_snow_veg)) DEALLOCATE (frac_snow_veg)
2766    IF ( ALLOCATED (frac_snow_nobio)) DEALLOCATE (frac_snow_nobio)
2767    IF ( ALLOCATED (snowrho)) DEALLOCATE (snowrho)
2768    IF ( ALLOCATED (snowgrain)) DEALLOCATE (snowgrain)
2769    IF ( ALLOCATED (snowtemp)) DEALLOCATE (snowtemp)
2770    IF ( ALLOCATED (snowdz)) DEALLOCATE (snowdz)
2771    IF ( ALLOCATED (snowheat)) DEALLOCATE (snowheat)
2772    IF ( ALLOCATED (cgrnd_snow)) DEALLOCATE (cgrnd_snow)
2773    IF ( ALLOCATED (dgrnd_snow)) DEALLOCATE (dgrnd_snow)
2774    IF ( ALLOCATED (lambda_snow)) DEALLOCATE(lambda_snow) 
2775    IF ( ALLOCATED (gtemp)) DEALLOCATE (gtemp)
2776    IF ( ALLOCATED (soilflxresid)) DEALLOCATE (soilflxresid)
2777    IF ( ALLOCATED (pgflux)) DEALLOCATE (pgflux)
2778    IF ( ALLOCATED (mc_layh)) DEALLOCATE (mc_layh)
2779    IF ( ALLOCATED (mcl_layh)) DEALLOCATE (mcl_layh)
2780    IF ( ALLOCATED (soilmoist)) DEALLOCATE (soilmoist)
2781    IF ( ALLOCATED (mc_layh_s)) DEALLOCATE (mc_layh_s)
2782    IF ( ALLOCATED (mcl_layh_s)) DEALLOCATE (mcl_layh_s)
2783!    IF ( ALLOCATED (mc_layh_pft)) DEALLOCATE (mc_layh_pft)
2784!    IF ( ALLOCATED (mcl_layh_pft)) DEALLOCATE (mcl_layh_pft)
2785!    IF ( ALLOCATED (soilmoist_pft)) DEALLOCATE (soilmoist_pft)
2786
2787!!pss:+
2788 !   IF ( ALLOCATED (fwet_out)) DEALLOCATE (fwet_out)
2789 !   IF ( ALLOCATED (drunoff_tot)) DEALLOCATE (drunoff_tot)
2790!!pss:-
2791
2792!! 3. Clear all allocated memory
2793
2794    CALL pft_parameters_clear
2795    CALL slowproc_clear 
2796    CALL diffuco_clear 
2797    CALL enerbil_clear 
2798    IF ( hydrol_cwrr ) THEN
2799       CALL hydrol_clear 
2800       CALL thermosoil_clear
2801    ELSE
2802       CALL hydrolc_clear 
2803       CALL thermosoilc_clear
2804    ENDIF
2805    CALL condveg_clear 
2806    CALL routing_clear
2807
2808  END SUBROUTINE sechiba_clear
2809
2810
2811!! ==============================================================================================================================\n
2812!! SUBROUTINE   : sechiba_var_init
2813!!
2814!>\BRIEF        Calculate air density as a function of air temperature and
2815!! pressure for each terrestrial pixel.
2816!!
2817!! RECENT CHANGE(S): None
2818!!
2819!! MAIN OUTPUT VARIABLE(S): air density (::rau, kg m^{-3}).
2820!!
2821!! REFERENCE(S) : None
2822!!
2823!! FLOWCHART    : None
2824!! \n
2825!_ ================================================================================================================================
2826
2827  SUBROUTINE sechiba_var_init (kjpindex, rau, pb, temp_air) 
2828
2829!! 0.1 Input variables
2830
2831    INTEGER(i_std), INTENT (in)                    :: kjpindex        !! Domain size - terrestrial pixels only (unitless)
2832    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: pb              !! Surface pressure (hPa)
2833    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: temp_air        !! Air temperature (K)
2834   
2835!! 0.2 Output variables
2836
2837    REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rau             !! Air density @tex $(kg m^{-3})$ @endtex
2838
2839!! 0.3 Modified variables
2840
2841!! 0.4 Local variables
2842
2843    INTEGER(i_std)                                 :: ji              !! Indices (unitless)
2844!_ ================================================================================================================================
2845   
2846!! 1. Calculate intial air density (::rau)
2847   
2848    DO ji = 1,kjpindex
2849       rau(ji) = pa_par_hpa * pb(ji) / (cte_molr*temp_air(ji))
2850    END DO
2851
2852    IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_var_init done '
2853
2854  END SUBROUTINE sechiba_var_init
2855
2856
2857!! ==============================================================================================================================\n
2858!! SUBROUTINE   : sechiba_end
2859!!
2860!>\BRIEF        Swap old for newly calculated soil temperature.
2861!!
2862!! RECENT CHANGE(S): None
2863!!
2864!! MAIN OUTPUT VARIABLE(S): soil temperature (::temp_sol; K)
2865!!
2866!! REFERENCE(S) : None
2867!!
2868!! FLOWCHART    : None
2869!! \n
2870!! ================================================================================================================================
2871
2872  SUBROUTINE sechiba_end (kjpindex, temp_sol_new, temp_sol_new_pft, temp_sol, temp_sol_pft)
2873                         
2874
2875!! 0.1 Input variables
2876
2877    INTEGER(i_std), INTENT (in)                       :: kjpindex           !! Domain size - terrestrial pixels only (unitless)
2878    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: temp_sol_new       !! New soil temperature (K)
2879    REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in)     :: temp_sol_new_pft       !! New soil temperature (K)
2880   
2881    !! 0.2 Output variables
2882    REAL(r_std),DIMENSION (kjpindex), INTENT (out)    :: temp_sol           !! Soil temperature (K)
2883    REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (out)    :: temp_sol_pft           !! Soil temperature (K)
2884
2885    !! 0.2 Local variables
2886    INTEGER(i_std)                                    :: ji,jv
2887!_ ================================================================================================================================
2888   
2889!! 1. Swap temperature
2890
2891    temp_sol(:) = temp_sol_new(:)
2892    DO jv = 1,nvm
2893        temp_sol_pft(:,jv) = temp_sol_new_pft(:,jv)
2894    ENDDO
2895
2896    IF (printlev_loc>=3) WRITE (numout,*) ' sechiba_end done '
2897
2898  END SUBROUTINE sechiba_end
2899
2900
2901
2902!! ==============================================================================================================================\n
2903!! SUBROUTINE   : sechiba_get_cmd
2904!!
2905!>\BRIEF       
2906!!
2907!! DESCRIPTION  :
2908!!
2909!! RECENT CHANGE(S):
2910!!
2911!! MAIN OUTPUT VARIABLE(S):
2912!!
2913!! REFERENCE(S) : None
2914!!
2915!! FLOWCHART    : None
2916!! \n
2917!! ================================================================================================================================
2918  SUBROUTINE sechiba_get_cmd(cmdin, src_rot, tgt_rot, prc_rot)
2919  ! 0.1 Input
2920  INTEGER(i_std),INTENT(in)  :: cmdin
2921  ! 0.2 Output
2922  INTEGER(i_std),INTENT(out) :: src_rot, tgt_rot
2923  REAL(r_std),INTENT(out)    :: prc_rot
2924
2925  ! 0.3 Local variables
2926
2927  !!!! --------------------------------------------------------------
2928    IF (cmdin > 1010000 .OR. cmdin < 0 ) THEN
2929        WRITE(numout,*) 'cmdin, ',cmdin
2930        STOP 'cmd error in sechiba_get_cmd'
2931    ENDIF
2932    IF (cmdin == 0) THEN
2933        tgt_rot = 0
2934        src_rot = 0
2935       
2936        prc_rot = 0.0
2937    ELSE
2938        tgt_rot = MOD(cmdin, 100)
2939        src_rot = MOD(FLOOR(FLOAT(cmdin)/100), 100)
2940
2941        prc_rot = FLOAT(FLOOR(FLOAT(cmdin)/10000))/100.0
2942        IF (printlev >=4) THEN
2943            WRITE(numout,*) 'xuhui: cmdin, tgt_rot, src_rot, prc_rot', cmdin, src_rot, tgt_rot, prc_rot
2944        ENDIF
2945        IF (prc_rot .GT. 1.0 .AND. prc_rot .LT. 1.0+0.01) THEN ! resolve potential  precision issues
2946            prc_rot = 1.0
2947        ENDIF
2948        !!! consistency check
2949        IF (prc_rot .GT. 1.0) THEN
2950            WRITE(numout,*) 'percent change larger than 1..., prc_rot',prc_rot
2951            STOP 'incorrect percent rotation, sechiba_get_cmd'
2952        ENDIF
2953        IF ( (tgt_rot .GT. nvm) .OR. ( .NOT. (ok_LAIdev(tgt_rot) .OR. (tgt_rot .EQ. 1)) ) ) THEN
2954            WRITE(numout,*) 'rotation target error: tgt_rot ', tgt_rot
2955            WRITE(numout,*) 'nvm, ok_LAIdev', nvm, ok_LAIdev
2956            STOP 'incorrect rotation target, sechiba_get_cmd'
2957        ENDIF
2958        IF ( (src_rot .GT. nvm) .OR. ( .NOT. (ok_LAIdev(src_rot) .OR. (src_rot .EQ. 1)) ) ) THEN
2959            WRITE(numout,*) 'rotation target error: src_rot ', src_rot
2960            WRITE(numout,*) 'nvm, ok_LAIdev', nvm, ok_LAIdev
2961            STOP 'incorrect rotation source, sechiba_get_cmd'
2962        ENDIF
2963    ENDIF
2964  END SUBROUTINE sechiba_get_cmd
2965
2966
2967!! ==============================================================================================================================\n
2968!! SUBROUTINE   : sechiba_interface_orchidee_inca
2969!!
2970!>\BRIEF        make the interface between surface and atmospheric chemistry
2971!!
2972!! DESCRIPTION  : This subroutine is called from INCA, the atmospheric chemistry model. It is used to transfer variables from ORCHIDEE to INCA.
2973!!
2974!! RECENT CHANGE(S): move from chemistry module to be more generic (feb - 2017)
2975!!
2976!! MAIN OUTPUT VARIABLE(S): emission COV to be transport by orchidee to inca in fields_out array
2977!!
2978!! REFERENCE(S) : None
2979!!
2980!! FLOWCHART    : None
2981!! \n
2982!! ================================================================================================================================
2983  SUBROUTINE sechiba_interface_orchidee_inca( &
2984       nvm_out, veget_max_out, veget_frac_out, lai_out, snow_out, &
2985       field_out_names, fields_out, field_in_names, fields_in)
2986
2987
2988    INTEGER, INTENT(out)                      :: nvm_out            !! Number of vegetation types
2989    REAL(r_std), DIMENSION (:,:), INTENT(out) :: veget_max_out      !! Max. fraction of vegetation type (LAI -> infty)
2990    REAL(r_std), DIMENSION (:,:), INTENT(out) :: veget_frac_out     !! Fraction of vegetation type (unitless, 0-1) 
2991    REAL(r_std), DIMENSION (:,:), INTENT(out) :: lai_out            !! Surface foliere
2992    REAL(r_std), DIMENSION (:)  , INTENT(out) :: snow_out           !! Snow mass [Kg/m^2]
2993
2994    !
2995    ! Optional arguments
2996    !
2997    ! Names and fields for emission variables : to be transport by Orchidee to Inca
2998    CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_out_names
2999    REAL(r_std),DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: fields_out
3000    !
3001    ! Names and fields for deposit variables : to be transport from chemistry model by INCA to ORCHIDEE.
3002    CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_in_names
3003    REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(IN)    :: fields_in
3004
3005
3006    ! Variables always transmitted from sechiba to inca
3007    nvm_out = nvm 
3008    veget_max_out(:,:)  = veget_max(:,:) 
3009    veget_frac_out(:,:) = veget(:,:) 
3010    lai_out(:,:)  = lai(:,:) 
3011    snow_out(:)  = snow(:) 
3012
3013    ! Call chemistry_flux_interface if at least one of variables field_out_names or
3014    ! field_in_names is present in the argument list of sechiba_interface_orchidee_inca when called from inca.
3015    IF (PRESENT(field_out_names) .AND. .NOT. PRESENT(field_in_names)) THEN
3016       CALL chemistry_flux_interface(field_out_names=field_out_names, fields_out=fields_out)
3017    ELSE IF (.NOT. PRESENT(field_out_names) .AND. PRESENT(field_in_names)) THEN
3018       CALL chemistry_flux_interface(field_in_names=field_in_names, fields_in=fields_in)
3019    ELSE IF (PRESENT(field_out_names) .AND. PRESENT(field_in_names)) THEN
3020       CALL chemistry_flux_interface(field_out_names=field_out_names, fields_out=fields_out, &
3021            field_in_names=field_in_names, fields_in=fields_in)
3022    ENDIF
3023
3024  END SUBROUTINE sechiba_interface_orchidee_inca
3025
3026
3027END MODULE sechiba
3028
Note: See TracBrowser for help on using the repository browser.