source: branches/publications/ORCHIDEE_GLUC_r6545/src_sechiba/sechiba.f90 @ 7346

Last change on this file since 7346 was 5089, checked in by albert.jornet, 6 years ago

Fix: missing OpenMP threadprivate declarations ( LMDZ+MICT )
Fix: group variable (stics module) was defined as real but passed as an integer
Fix: masec variable (stics module) it was modified but passed as input argument only

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