source: branches/ORCHIDEE_3_CMIP6/ORCHIDEE/src_sechiba/sechiba.f90 @ 7599

Last change on this file since 7599 was 7245, checked in by nicolas.vuichard, 3 years ago

improve Carbon mass balance closure. See ticket #785

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