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

Last change on this file since 8566 was 8557, checked in by josefine.ghattas, 8 weeks ago

Added print for pixel and lat/lon as done in the trunk. This is usefull information printed only once.

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