source: branches/publications/ORCHIDEE_Biochar/src_sechiba/sechiba.f90 @ 7746

Last change on this file since 7746 was 7366, checked in by simon.bowring, 3 years ago

Biochar version

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