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

Last change on this file since 7525 was 7525, checked in by sebastiaan.luyssaert, 2 years ago

Fix for ticket #816

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