source: branches/publications/ORCHILEAK-Gommet/src_sechiba/sechiba.f90 @ 7746

Last change on this file since 7746 was 5315, checked in by ronny.lauerwald, 7 years ago

Bug fix DOC inputs in floodplains and swamps

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