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

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

Added routing_simple module and routing_wrapper. No changes in defalut set up.

See ticket #581

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