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

Last change on this file since 7326 was 7326, checked in by josefine.ghattas, 3 years ago

Corrected bug on carbon balance closure. See ticket #785
Integration in branch 2_2 done by P. Cadule

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