source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_sechiba/sechiba.f90 @ 7541

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