source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_sechiba/sechiba.f90 @ 8398

Last change on this file since 8398 was 7656, checked in by christophe.dumas, 2 years ago

Delete the unused subroutine thermosoi_ice and the unused transmitted variables in thermosoil

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