source: branches/publications/ORCHIDEE-PEAT_r5488/src_sechiba/sechiba.f90 @ 8787

Last change on this file since 8787 was 5488, checked in by chunjing.qiu, 6 years ago

C balance checked

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