Context Navigation

hydraulic_arch.f90 @ 7852

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

Enhanced consistency of variable names: input has been changed in n_input throughout the code and the variable name vegstress introduced in sechiba is now also used in stomate. Enhnaced computational consistency: Pgap_cumul is used in stomate rather than recalculating it before calculating light_tran_to_floor_season. Edited getin_p while checking the code (but no real changes were made) and added several missing stomate and sechiba variables to age_class_distr to improve the 1+1=2 issue when comparing a model run with against a run without age classes. Finally: Write warning 10b in allocation to the history file rather than the out_orchidee file

File size: 126.5 KB

Line
1	!==================================================================================================================================
2	! MODULE : hydraulic_arch
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF This module calculates the amount of water that is available for the leaves for transpiration,
10	! accounting for effects of hydraulic architecture following (Hickler et al. 2006)
11	!!
12	!!\n DESCRIPTION : (1) Calculate soil water potential from soil water content as calculated in hydrol.f90
13	!! according to van Genuchten (1980). Soil water potential is weighted by the relative share of the
14	!! of theroots in each layer.
15	!! (2) Calculate pressure difference between leaves and soil (Whitehead 1998)
16	!! (3) Calculate root (Weatherly 1998), sapwood (Magnani et al. 2000, Tyree et al. 1994,
17	!! Sperry et al. 1998) and leave resistances (Sitch et al. 2003.
18	!! (4) Calculate the supply of water for transpiration by the leaves using Darcy's law (Whitehead 1998)
19	!!
20	!! RECENT CHANGE(S):(1) Soil to root resistance has been added to the calculations of the soil water potential weighted
21	!! by root density. The development was made by Emilie Joetzjer based on
22	!! Williams et al., 2001 (tree physiology) and Duursma et al., 2012 (GMD). Also see Bonan et al. 2014 (GMD)
23	!!
24	!! REFERENCE(S) : None
25	!!
26	!! SVN :
27	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_sechiba/sechiba.f90 $
28	!! $Date: 2013-05-06 14:14:36 +0200 (Mon, 06 May 2013) $
29	!! $Revision: 1295 $
30	!! \n
31	!_ ================================================================================================================================
32
33	MODULE hydraulic_arch
34
35	USE ioipsl
36	USE xios_orchidee
37
38	! modules used :
39	USE constantes
40	USE constantes_soil
41	USE constantes_soil_var
42	USE pft_parameters
43	USE ioipsl_para
44	USE function_library, ONLY: wood_to_height_eff, cc_to_lai, &
45	biomass_to_coupled_lai, biomass_to_lai, &
46	Arrhenius, Arrhenius_modified, &
47	Arrhenius_modified_nodim
48	USE qsat_moisture
49	USE mleb
50	USE enerbil, ONLY: enerbil_main, enerbil_write
51
52	IMPLICIT NONE
53
54	! Private and public routines
55	PRIVATE
56	PUBLIC hydraulic_arch_main, mleb_diffuco_trans_co2_short
57
58
59	LOGICAL, SAVE :: l_first_sechiba = .TRUE. !! Flag controlling the intialisation (true/false)
60	!$OMP THREADPRIVATE(l_first_sechiba)
61	INTEGER(i_std), SAVE :: printlev_loc !! Local printlev in slowproc module
62	!$OMP THREADPRIVATE(printlev_loc)
63	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: nvan !! VanGenuchten coeficients n (unitless) ! RK:
64	!$OMP THREADPRIVATE(nvan)
65	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: avan !! Van Genuchten coeficients a!@tex!$(mm^{-1})$!@endtex
66	!$OMP THREADPRIVATE(avan)
67	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcr !! Residual volumetric water content !! !@tex !$(m^{3} !m^{-3})$!@endtex
68	!$OMP THREADPRIVATE(mcr)
69	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcs !! Saturated volumetric water content !@tex !$(m^{3} !m^{-3})$!@endtex
70	!$OMP THREADPRIVATE(mcs)
71
72	CONTAINS
73
74
75
76	!! =============================================================================================================================
77	!! SUBROUTINE: hydraulic_arch_main
78	!!
79	!>\BRIEF: call the relevant modules for the multi-layer energy budget and hydraulic stress
80	!!
81	!! DESCRIPTION: call the relevant modules for the multi-layer energy budget and hydraulic stress
82	!! This subroutine was preivously stored in sechiba.f90 and was named sechiba_mleb_hs
83	!! \n
84	!_ ==============================================================================================================================
85	SUBROUTINE hydraulic_arch_main( &
86	kjit, kjpij, kjpindex, rest_id, &
87	rest_id_stom, hist_id, hist2_id, hist_id_stom, &
88	hist_id_stom_IPCC, index, indexveg, ldrestart_read, &
89	ldrestart_write, njsc, u, v, &
90	qair, precip_rain, precip_snow, lwdown, &
91	swnet, swdown, temp_air, &
92	petAcoef, peqAcoef, petBcoef, peqBcoef, &
93	pb, tq_cdrag, epot_air, sinang, &
94	ccanopy, emis, soilcap, soilflx, &
95	rau, temp_growth, vbeta23, humrel, &
96	lai, max_height_store, gs_distribution, gs_diffuco_output, &
97	gstot_component, gstot_frac, assim_param, circ_class_biomass, &
98	circ_class_n, ksoil, lai_per_level, laieff_isotrop, &
99	Light_Abs_Tot, Light_Tran_Tot, qsintmax, qsintveg, &
100	snowdz, stempdiag, swc, &
101	vbeta1, vbeta2, vbeta4, vbeta5, &
102	veget_max, veget, z_array_out, &
103	temp_sol_new, tsol_rad, temp_sol_add, vevapflo, &
104	vevapnu, vevapsno, stressed, unstressed, &
105	transpot, vevapwet, transpir, &
106	e_frac, &
107	vbeta2sum, vbeta3sum, fluxsens, fluxlat, &
108	vevapp, vbeta, evapot, evapot_corr, &
109	qsurf, temp_sol, pgflux, u_speed, &
110	vbeta3, vbeta3pot, rveget, rstruct, &
111	cimean, gsmean, gpp, transpir_supply, &
112	leaf_ci, warnings, JJ_out, assimi_lev, &
113	cresist, info_limitphoto, vessel_loss, root_profile)
114
115	!! 0.1 Input variables
116	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
117	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
118	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
119	INTEGER(i_std), INTENT(in) :: rest_id !! _Restart_ file identifier (unitless)
120	INTEGER(i_std), INTENT(in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
121	INTEGER(i_std), INTENT(in) :: hist_id !! _History_ file identifier (unitless)
122	INTEGER(i_std), INTENT(in) :: hist2_id !! _History_ file 2 identifier (unitless)
123	INTEGER(i_std), INTENT(in) :: hist_id_stom !! STOMATE's _History_ file identifier
124	INTEGER(i_std), INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file identifier
125	INTEGER(i_std), DIMENSION(kjpindex), INTENT (in) :: index !! Indices of the pixels on the map.
126	INTEGER(i_std), DIMENSION(kjpindex*nvm), INTENT(in) :: indexveg !! Indeces of the points on the 3D map
127	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
128	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
129	INTEGER(i_std), DIMENSION(kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class
130	!! in the grid cell used for hydrology, (1-nscm, unitless)
131	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: u !! Lowest level wind speed in direction u
132	!! @tex $(m.s^{-1})$ @endtex
133	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: v !! Lowest level wind speed in direction v
134	!! @tex $(m.s^{-1})$ @endtex
135	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
136	!! @tex $(kg kg^{-1})$ @endtex
137	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
138	!! @tex $(kg m^{-2})$ @endtex
139	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
140	!! @tex $(kg m^{-2})$ @endtex
141	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux
142	!! @tex $(W m^{-2})$ @endtex
143	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux
144	!! @tex $(W m^{-2})$ @endtex
145	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux
146	!! @tex $(W m^{-2})$ @endtex
147	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
148	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: petAcoef !! Coefficients A for T from the Planetary
149	!! Boundary Layer
150	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: peqAcoef !! Coefficients A for q from the Planetary
151	!! Boundary Layer
152	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: petBcoef !! Coefficients B for T from the Planetary
153	!! Boundary Layer
154	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: peqBcoef !! Coefficients B for q from the Planetary
155	!! Boundary Layer
156	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
157	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: tq_cdrag !! Surface drag coefficient (-)
158	!! @tex $(m.s^{-1})$ @endtex
159	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: epot_air !! Air potential energy (??J)
160	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: sinang !! Sine of the solar angle (unitless)
161	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: ccanopy !! CO2 concentration in the canopy (ppm)
162	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: emis !!
163	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity (J K^{-1])
164	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: soilflx !! Soil flux (W m^{-2})
165	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3})
166	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
167	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in) :: vbeta23 !! Beta for fraction of wetted foliage that will
168	!! transpire once intercepted water has evaporated (-)
169	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture availability [0-1]
170	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !!
171	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: max_height_store !!
172	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in) ::assimi_lev
173	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot+1), INTENT(in)::info_limitphoto !! Stored for mleb_diffuco_
174	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in) ::JJ_out
175	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in) :: gs_distribution !!
176	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in) :: gs_diffuco_output !!
177	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in) :: gstot_component !!
178	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in) :: gstot_frac !!
179	REAL(r_std), DIMENSION(kjpindex,nvm,npco2), INTENT(in) :: assim_param !!
180	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc,nparts,nelements), INTENT(in) :: circ_class_biomass
181	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc), INTENT(in) :: circ_class_n
182	REAL(r_std), DIMENSION(kjpindex,nslm,nstm), INTENT (in) :: ksoil !! Hydraulic conductivity mm/d
183	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in) :: lai_per_level !! This is the LAI per vertical level
184	REAL(r_std), DIMENSION(kjpindex,nlevels_tot,nvm), INTENT(in) :: laieff_isotrop !! Effective LAI
185	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT (in) :: Light_Abs_Tot !! Absorbed radiation per level for photosynthesis
186	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT (in) :: Light_Tran_Tot !! Absorbed radiation per level for photosynthesis
187	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation
188	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
189	REAL(r_std), DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer [m]
190	REAL(r_std), DIMENSION(kjpindex,nslm),INTENT(in) :: stempdiag !! Soil temperature (K)
191	REAL(r_std), DIMENSION(kjpindex,nslm,nstm), INTENT (in) :: swc !! Volumetric soil water content
192	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-)
193	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in) :: vbeta2 !! Interception resistance (-)
194	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: vbeta4 !! Bare soil resistance (-)
195	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance
196	REAL(r_std), DIMENSION(kjpindex, nvm), INTENT(in) :: veget_max !! Max. fraction of vegetation type (LAI -> infty, unitless)
197	REAL(r_std), DIMENSION(kjpindex, nvm), INTENT(in) :: veget !! Vegetation fraction for each pft
198	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc,nlevels_tot), INTENT(in) :: z_array_out !! An output of h_array, to use in sechiba
199	REAL(r_std),DIMENSION (kjpindex,nvm,nslm), INTENT(in) :: root_profile !! Normalized root mass/length fraction in each soil layer
200	!! (0-1, unitless)
201
202	!! 0.2 Output variables
203	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: cresist !! coefficient for resistances (??)
204	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: temp_sol_new !! New ground temperature (K)
205	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: tsol_rad !! Radiative surface temperature (W/m2)
206	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
207	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: vevapflo !! Floodplains evaporation
208	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: vevapnu !! Bare soil evaporation (mm day^{-1})
209	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: vevapsno !! Snow evaporation (mm day^{-1})
210	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: stressed !! GPP when the PFT is stressed (g C m-2 dt-1)
211	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: unstressed !! GPP when the PFT is unstressed (g C m-2 dt-1)
212	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: transpot !! Potential transpiration
213	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: vevapwet !! Interception (mm day^{-1})
214	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: transpir !! Transpiration (mm day^{-1})
215	REAL(r_std), DIMENSION(kjpindex,nvm,nslm,nstm), INTENT(out) :: e_frac !! Fraction of water transpired supplied by individual layers (no units)
216	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: vessel_loss !! Proportion of conductivity lost due to cavitation in the xylem (no unit)
217
218	!! 0.3 Modified variables
219	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: vbeta2sum !! sum of vbeta2 coefficients across all PFTs (-)
220	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: vbeta3sum !! sum of vbeta3 coefficients across all PFTs (-)
221	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: fluxsens !! Sensible heat flux
222	!! @tex $(W m^{-2})$ @endtex
223	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: fluxlat !! Latent heat flux
224	!! @tex $(W m^{-2})$ @endtex
225	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: vevapp !! Total of evaporation
226	!! @tex $(kg m^{-2} days^{-1})$ @endtex
227	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: vbeta !! Resistance coefficient (-)
228	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: evapot !!
229	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: evapot_corr !!
230	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: qsurf !!
231	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: temp_sol !!
232	REAL(r_std), DIMENSION(kjpindex), INTENT(inout) :: pgflux !! Net energy into snowpack(W/m^2)
233	REAL(r_std), DIMENSION(jnlvls), INTENT(inout) :: u_speed !! Canopy wind speed profile - see comment in mleb_finalize
234	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: vbeta3 !!
235	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: vbeta3pot !!
236	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: rveget !!
237	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: rstruct !!
238	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: cimean !!
239	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: gsmean !!
240	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: gpp !! gC m^{-2}pixel
241	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: transpir_supply !! Supply of water for transpiration (mm timestep^{-1})
242	REAL(r_std), DIMENSION(kjpindex,nvm,nlai), INTENT(inout) :: leaf_ci !! intercellular CO2 concentration (ppm)
243	REAL(r_std), DIMENSION(kjpindex,nvm,nwarns), INTENT(inout) :: warnings !! A warning counter
244
245
246	!! 0.4 Local variables
247	! variables introduced as a consequence of running enerbil twice (if needed, as a result of restriction of transpiration by
248	! the supply term). If running for the second time, we should start with the same input variables as the first run
249
250	INTEGER(i_std) :: ji, jv
251	REAL(r_std), DIMENSION (kjpindex) :: evapot_save !! Soil Potential Evaporation (mm day^{-1})
252	REAL(r_std), DIMENSION (kjpindex) :: evapot_corr_save !! Soil Potential Evaporation Correction (mm day^{-1})
253	REAL(r_std), DIMENSION (kjpindex) :: diffevap !! Difference betwence vevapp and composing fluxes (Kg/m^2/s)
254	REAL(r_std), DIMENSION (kjpindex) :: temp_sol_save !! Soil temperature (K)
255	REAL(r_std), DIMENSION (kjpindex) :: qsurf_save !! Surface specific humidity (kg kg^{-1})
256	REAL(r_std), DIMENSION (kjpindex) :: tsol_rad_save !! Tsol_rad (W m^{-2})
257	REAL(r_std), DIMENSION (kjpindex) :: netrad_save !! Net radiation flux (W m^{-2})
258	REAL(r_std), DIMENSION (kjpindex) :: vbeta_save !! Resistance coefficient (-)
259	REAL(r_std), DIMENSION (kjpindex) :: vbeta2sum_save !! sum of vbeta2 coefficients across all PFTs (-)
260	REAL(r_std), DIMENSION (kjpindex) :: vbeta3sum_save !! sum of vbeta3 coefficients across all PFTs (-)
261	REAL(r_std), DIMENSION (kjpindex,nvm) :: vbeta3_save !! Beta for Transpiration (unitless)
262	REAL(r_std), DIMENSION (kjpindex,nvm) :: vbeta3pot_save !! Beta for Potential Transpiration
263	REAL(r_std), DIMENSION (kjpindex) :: pgflux_save !! Net energy into snowpack(W/m^2)
264	REAL(r_std), DIMENSION (kjpindex) :: qsol_sat_new !! New saturated surface air moisture (kg kg^{-1})
265	REAL(r_std), DIMENSION (kjpindex) :: qair_new
266	REAL(r_std), DIMENSION (kjpindex) :: vbeta3_save_test
267	REAL(r_std), DIMENSION (kjpindex) :: netrad !! Net radiation (W m^{-2})
268	REAL(r_std), DIMENSION (kjpindex) :: lwabs !! LW radiation absorbed by the surface (W m^{-2})
269	REAL(r_std), DIMENSION (kjpindex) :: lwnet !! Net Long-wave radiation (W m^{-2})
270	REAL(r_std), DIMENSION (kjpindex) :: fluxsubli !! Energy of sublimation (mm day^{-1})
271
272	REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_vbeta3
273	REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_rveget
274	REAL(r_std), DIMENSION (nlevels_tot,kjpindex,nvm) :: transpir_supply_column !! (SECHIBA TOP LEVEL) Supply of water for transpiration
275	!+++CHECK+++
276	! Variable dimensions are for a single pixel. Causing 1+1 problems
277	! when running a larger domain. Needs to be corrected when implementing
278	! a global use of the multi-layer energy budget.
279	!!$ REAL(r_std), DIMENSION (nlevels_tot) :: Light_Abs_Tot_mean !! (SECHIBA TOP LEVEL) total light absorption for a given canopy level
280	!!$ REAL(r_std), DIMENSION (nlevels_tot) :: Light_Alb_Tot_mean !! (SECHIBA TOP LEVEL) total albedo for a given level
281	!+++++++++++
282
283	LOGICAL :: hydrol_flag !! flag that 'trips' the energy budget for each grid square
284	LOGICAL, DIMENSION (kjpindex) :: hydrol_flag2 !! flag that 'trips' the energy budget for each grid square
285	LOGICAL, DIMENSION (nbp_glo) :: hydrol_flag2_glo !! hydrol_flag2 on the full global domain for all processors
286	LOGICAL, DIMENSION (kjpindex,nvm) :: hydrol_flag3 !! flag that 'trips' the energy budget for each grid square and PFT
287
288	LOGICAL, DIMENSION (kjpindex,nlevels_tot) :: hydrol_flag2_column !! flag that 'trips' the energy budget for each grid square, and canopy level
289	LOGICAL, DIMENSION (kjpindex,nvm,nlevels_tot) :: hydrol_flag3_column !! flag that 'trips' the energy budget for each grid square, canopy level and PFT
290
291	REAL(r_std), DIMENSION (kjpindex,nvm) :: ratio_vbeta3 !! Reduction rate of vbeta3 when transpir_supply
292	!! is smaller than demand
293
294	REAL(r_std), DIMENSION(kjpindex,nvm) :: cresist_save !!temporal
295
296	! ==============================================================================================================================
297
298	! Define a set of flags, required to signal if we need to recalculate the
299	! energy budget depending on the hydraulic stress being a limiting factor to the
300	! latent heat flux.
301	hydrol_flag = .FALSE. ! trip flag per cycle (no dimension)
302	hydrol_flag2(:) = .FALSE. ! trip flag per grid square
303	hydrol_flag3(:,:) = .FALSE. ! trip flag per grid square and PFT
304
305	! Store vbeta3 for the calculation of ratio_vbeta3
306	vbeta3_save(:,:) = vbeta3(:,:)
307	ratio_vbeta3(:,:) = un
308
309	! Store values of inout variables before calculating the energy budget
310	! so it can be recalculated with the exact same initial conditions when
311	! needed.
312	evapot_save(:) = evapot(:)
313	evapot_corr_save(:) = evapot_corr(:)
314	temp_sol_save(:) = temp_sol(:)
315	qsurf_save(:) = qsurf(:)
316	tsol_rad_save(:) = tsol_rad(:)
317	pgflux_save(:) = pgflux(:)
318
319
320	IF (.not.ok_mleb) THEN
321
322	! Calculate the energy budget for the first time. The main reason we need
323	! this calculation is to have an estimate of transpir.
324	CALL enerbil_main (kjit, kjpindex, netrad, lwabs, lwnet, fluxsubli, &
325	& index, indexveg, lwdown, swnet, &
326	& epot_air, temp_air, u, v, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, &
327	& pb, rau, vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
328	& emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
329	& vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo,&
330	& temp_sol, tsol_rad,temp_sol_new, qsurf, evapot, evapot_corr, &
331	& precip_rain, pgflux, snowdz, temp_sol_add, qair_new, qsol_sat_new)
332
333	ELSE
334
335	!+++CHECK+++
336	! Add the inout variables for mleb
337	CALL ipslerr_p(3,'hydraulic_arch_main','The recalculation of the energybudget',&
338	'needs to be checked when mleb rather than enerbil_main is used',&
339	'Save the intent(inout) variables for mleb')
340	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
341	! when running a larger domain. Needs to be corrected when implementing
342	! a global use of the multi-layer energy budget.
343	!+++++++++++
344
345	CALL mleb_main (kjit, kjpindex, ldrestart_read, ldrestart_write, &
346	& index, indexveg, lwdown, swnet, swdown, epot_air, temp_air, &
347	& u, v, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, pb, rau, &
348	& vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
349	& emis, soilflx, soilcap, tq_cdrag, humrel, &
350	& fluxsens, fluxlat, vevapp, transpir, transpot, vevapnu, vevapwet, &
351	& vevapsno, vevapflo, temp_sol, tsol_rad, temp_sol_new, qsurf, &
352	& evapot, evapot_corr, rest_id, hist_id, hist2_id, &
353	& ok_mleb_history_file, &
354	& transpir_supply, hydrol_flag, hydrol_flag2, hydrol_flag3, vbeta2sum, &
355	& vbeta3sum, veget_max, qsol_sat_new, qair_new, &
356	& veget, &
357	!!$ & Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
358	& laieff_isotrop, &
359	& z_array_out, transpir_supply_column, u_speed, &
360	& profile_vbeta3, profile_rveget, max_height_store, &
361	& precip_rain, pgflux, snowdz, temp_sol_add)
362
363	END IF
364
365	! Store the proxy for unstressed ecosystem functioning. This could be gpp, transpiration,
366	! evaporation. They all differ conceptually and numerically but the idea is the same.
367	! In an initial approach we tried to calculate the stress as a ratio of the initial
368	! transpiration over water supplied through the plant. This approach resulted in
369	! inconsistencies likely because the many feedbacks in the calculation of transpir. For
370	! this reason it was decided to drop this approach and use the ratio between the initial
371	! and adjusted gpp instead. This ratio will be used in stomate to adjust the allocation.
372	! Note that the feedback on gpp is calculated based on the ratio between ::transpir_supply
373	! and transpir at the half hourly time step. In stomate we use the ratio between stressed
374	! and unstressed but at the daily time step. Night time values should be ignored in the
375	! daily calculation. If night time values are considered to unstressed this will skew
376	! the result.
377	unstressed(:,:) = gpp(:,:)
378
379	IF (ok_hydrol_arch) THEN
380
381	CALL hydraulic_arch_calc (kjpindex, circ_class_biomass, circ_class_n, transpir, stempdiag, &
382	temp_air, swc, transpir_supply, njsc, veget, &
383	transpir_supply_column, ksoil, e_frac, vessel_loss, root_profile)
384
385	DO ji = 1, kjpindex
386
387	!don't restrict gs during night
388	DO jv = 1, nvm
389	IF (sinang(ji) .LT. min_sechiba) THEN
390	! Night time we will never recalculate the energy budget
391	hydrol_flag2(ji) = .FALSE.
392	hydrol_flag3(ji,jv) = .FALSE.
393	ELSE
394	! The ratio between ::transpir_supply and ::transpir is used to determine
395	! whether we will recalculate gpp, transpir, ... (thus enerbil). This
396	! is done every half hour. The waterstress that will be passed to stomate
397	! to adjust the allocation will be calculated at the daily time step.
398	! It is important NOT to include night time values in the daily mean
399	! because more southern sides have shorter days during the growing season
400	! the number of half hours with day light biases the waterstress if it
401	! is assumed that there is no waterstress at night.
402	! Note that if transpir_suppy = zero this is likely caused by a soil for
403	! which all 11 layers are below wilting point. If this is the case it is
404	! correct to have transpir_supply = 0 and thus gpp = 0. This behaviour is
405	! not frequently observed. Even not in very dry regions. It suggests that
406	! the plant response to soil drying is too slow, implying transpir_supply
407	! is too high the days/week before the drought. Indeed the scheme
408	! to calculate the maximum transpir_supply is static, i.e. it always calculates
409	! the maximum transpir_supply (and therefore stomata might close too late or too
410	! little).
411	IF (veget_max(ji,jv) .GT. zero .AND. &
412	veget(ji,jv) .GT. zero) THEN
413	! comparing the transpiration in tree-stand scale
414	IF ( transpir_supply(ji,jv) .LT. transpir(ji,jv) ) THEN
415	IF(printlev_loc>=3) WRITE(numout,*) 'hydrol_flag2&3 set to TRUE, pft is: ', jv
416	hydrol_flag2(ji) = .TRUE.
417	hydrol_flag3(ji,jv) = .TRUE.
418	ELSE
419	IF(printlev_loc>=3) WRITE(numout,*) 'hydrol_flag2 set to FALSE, pft is: ', jv
420	END IF
421	ELSE
422	! There is no vegetation here, so we can skip it. The flag at the
423	! pixel level (::hydrol_flag2) should not change its value.
424	hydrol_flag3(ji,jv) = .FALSE.
425	ENDIF
426	END IF ! sinang(ji) .LT. min_sechiba
427	END DO ! jv = 1, nvm
428	END DO ! ji = 1, kjpindex
429
430	! The energy budget is only calculated again, if the above conditions are satisfied
431	! for any of the grid squares, and the feedback of waterstress on stomatal closure
432	! (ok_gs_feedback) is activated
433	IF (ok_gs_feedback) THEN
434	! Gather hydrol_flag2 on the full global domain and set hydrol_flag true
435	! if at least one grid cell contains hydrol_flag2()=true.
436	CALL gather(hydrol_flag2, hydrol_flag2_glo)
437	IF (is_root_prc) THEN
438	hydrol_flag=.FALSE.
439	DO ji = 1, nbp_glo
440	IF (hydrol_flag2_glo(ji)) THEN
441	hydrol_flag = .TRUE.
442	EXIT
443	END IF
444	END DO
445	END IF
446	! Copy hydrol_flag to all processors
447	CALL bcast(hydrol_flag)
448	ELSE
449	! If the feedback of water stress on stomatal closure is not activated,
450	! then hydrol_flag is set to false, and no energy re-calculations are conducted
451	hydrol_flag=.FALSE.
452	END IF
453
454	IF (hydrol_flag) THEN
455	DO ji = 1, kjpindex
456	IF (hydrol_flag2(ji)) THEN
457	DO jv = 1, nvm
458	IF (.NOT. hydrol_flag3(ji, jv) .OR. lai(ji,jv) .LT. 0.5 ) THEN
459	! this is for cases when the supply term is higher than the transpiration
460	! i.e. the normal state of affairs. The condition for lai avoids that
461	! we are getting crazy values when recalculating the gpp. Towards the
462	! end of the growing season the lai is decreasing which results in a very
463	! high R_shoot. If there is serious water stress at the time that the lai
464	! is very low, the recalculated gpp values go through the roof. We use
465	! this arbitrary threshold to avoid this behaviour. This implies that small
466	! canopies will not experience water stress in the model.
467	ELSE
468	! To calculate the gs_effective (the effective stomatal conductance), and
469	! GPP first we need to calculate a constrained vbeta3, and work back from
470	! there, using a tweaked version of diffuco_trans_co2. We calculate the
471	! vbeta3 based on the ::transpir_supply term for each grid and pft. We should
472	! keep the spatial weighting fraction,::veget_max, for vebta3. Even if the
473	! calculation of ::transpir_supply is at the stand scale.
474	IF (veget_max(ji,jv) .GT. zero .AND. veget(ji,jv) .GT. zero) THEN
475	vbeta3(ji,jv) = transpir_supply(ji,jv) / &
476	& ( dt_sechiba * (un - vbeta1(ji)) * (qsol_sat_new(ji) - qair_new(ji)) * &
477	& ( rau(ji) * (MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))) * tq_cdrag(ji)))
478	ELSE
479	! No vegetation present for this condition, so vbeta3 is set equal to zero.
480	vbeta3(ji,jv) = zero
481	ENDIF
482
483	! Calculate the ratio between the stressed and unstressed vbeta3. It will be
484	! used to correct some extreme values.
485	ratio_vbeta3(ji,jv) = vbeta3(ji,jv)/vbeta3_save(ji,jv)
486
487	! Recalculate leaf_ci, gsmean and gpp for the new vbeta3 that matches the
488	! water limitation
489	CALL mleb_diffuco_trans_co2_short(assimi_lev, assim_param, ccanopy, gstot_frac, ji, jv, &
490	kjpindex, Light_Abs_Tot, Light_Tran_Tot, lai, lai_per_level, pb, qair, tq_cdrag, &
491	qsintmax, qsintveg, ratio_vbeta3, rstruct, rveget, swdown, temp_air, temp_growth, &
492	u, v, vbeta23, vbeta3, veget_max, veget, info_limitphoto, JJ_out, cimean, cresist, &
493	leaf_ci, vbeta3pot, gsmean, gpp)
494
495	! Clean some unexpected/extreme values (less than 1% of the cases)
496	IF(gpp(ji,jv) .EQ. unstressed(ji,jv)) THEN
497	gpp(ji,jv) = ratio_vbeta3(ji,jv) * unstressed(ji,jv)
498	ELSE IF( gpp(ji,jv) .LE. zero) THEN
499	gpp(ji,jv) = zero
500	ELSE IF( gpp(ji,jv) .LE. 0.1ratio_vbeta3(ji,jv)unstressed(ji,jv)) THEN
501	gpp(ji,jv) = 0.1ratio_vbeta3(ji,jv)unstressed(ji,jv)
502	ELSE IF( gpp(ji,jv) .GT. 100 ) THEN
503	IF (err_act.GT.1) THEN
504	WRITE(numout,*)'GPP',gpp(ji,jv)
505	WRITE(numout,*)'vbeta3',vbeta3(ji,jv)
506	WRITE(numout,*)'vbeta3_save',vbeta3_save(ji,jv)
507	WRITE(numout,*)'transpir_supply',transpir_supply(ji,jv)
508	WRITE(numout,*)'vbeta1',vbeta1(ji)
509	WRITE(numout,*)'qsol_sat_new',qsol_sat_new(ji)
510	WRITE(numout,*)'qair_new',qair_new(ji)
511	WRITE(numout,*)'rau',rau(ji)
512	WRITE(numout,*)'tq_cdrag',tq_cdrag(ji)
513	WRITE(numout,*)'transpir',transpir(ji,jv)
514	WRITE(numout,*)'rstruct',rstruct(ji,jv)
515	WRITE(numout,*)'rveget',rveget(ji,jv)
516	WRITE(numout,*)'gsmean',gsmean(ji,jv)
517	CALL ipslerr_p(2,'hydraulic_arch_main', &
518	'The recalculation of GPP is too high(>100)!',&
519	'If this happens a lot, you should take a closer look','')
520	ENDIF
521	ENDIF
522	END IF ! (hydrol_flag3(ji, jv) .eq. .FALSE.)
523	ENDDO ! jv = 1, nvm
524	ELSE
525	! do nothing
526	END IF ! (hydrol_flag2(ji) .eq. .TRUE.)
527	END DO ! ji = 1, kjpindex
528
529
530	! mleb_main and enerbill_main have to be called here again
531	! (to recalculate the energy budget because of the
532	! supply term restriction.
533
534	! Before 2nd call to enerbil, recall initial values from before
535	! mleb was called the first time.
536	evapot(:) = evapot_save(:)
537	evapot_corr(:) = evapot_corr_save(:)
538	temp_sol(:) = temp_sol_save(:)
539	qsurf(:) = qsurf_save(:)
540	tsol_rad(:) = tsol_rad_save(:)
541	pgflux(:) = pgflux_save(:)
542
543	IF (.not. ok_mleb) THEN
544
545	! Second call calculates the energy budget with the water limited value of vbeta3
546	CALL enerbil_main (kjit, kjpindex, netrad, lwabs, lwnet, &
547	fluxsubli, index, indexveg, lwdown, swnet, &
548	epot_air, temp_air, u, v, petAcoef, petBcoef, qair, peqAcoef,peqBcoef,&
549	pb, rau, vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
550	emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
551	vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo,&
552	temp_sol, tsol_rad,temp_sol_new, qsurf, evapot, evapot_corr, precip_rain, &
553	pgflux, snowdz, temp_sol_add, qair_new, qsol_sat_new)
554
555	ELSE
556
557	!+++CHECK+++
558	! (1) Add the inout variables for mleb
559	! (2) Variable dimensions of xxx_Tot_mean are for a single pixel. Causing 1+1 problems
560	! when running a larger domain. Needs to be corrected when implementing
561	! a global use of the multi-layer energy budget.
562	!+++++++++++
563	CALL mleb_main (kjit, kjpindex, ldrestart_read, ldrestart_write, &
564	& index, indexveg, lwdown, swnet, swdown, epot_air, temp_air, &
565	& u, v, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, pb, rau, &
566	& vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
567	& emis, soilflx, soilcap, tq_cdrag, humrel, &
568	& fluxsens, fluxlat, vevapp, transpir, transpot, vevapnu, vevapwet, &
569	& vevapsno, vevapflo, temp_sol, tsol_rad, temp_sol_new, qsurf, &
570	& evapot, evapot_corr, rest_id, hist_id, hist2_id, &
571	& ok_mleb_history_file, &
572	& transpir_supply, hydrol_flag, hydrol_flag2, &
573	& hydrol_flag3, vbeta2sum, vbeta3sum, veget_max, qsol_sat_new, qair_new, &
574	& veget, &
575	!!$ & Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
576	& laieff_isotrop, z_array_out, transpir_supply_column, &
577	& u_speed, profile_vbeta3, profile_rveget, max_height_store, &
578	& precip_rain, pgflux, snowdz, temp_sol_add)
579
580	END IF
581
582	ELSE
583	! do nothing
584	END IF ! (hydrol_flag .eq. .TRUE.)
585
586	END IF ! IF (ok_hydrol_arch) THEN
587
588	! Make sure there is always a stressed and an unstressed value to avoid numerical
589	! issues in stomate. Store the proxy for stressed ecosystem functioning. This value
590	! is used in stomate to calculate the waterstress used in allocation. See comment
591	! where unstressed(:,:) is defined.
592	stressed(:,:) = gpp(:,:)
593
594
595	!! 2.3.1 Write enerbil data to the history files.
596	IF (ok_mleb) THEN
597	CALL mleb_write (kjit, kjpindex, index, lwdown, temp_sol, &
598	temp_sol_new, evapot, evapot_corr, hist_id, hist2_id, &
599	vevapp, vevapwet, transpir, vevapnu, vevapsno, vevapflo)
600	ELSE
601	CALL enerbil_write (kjit, kjpindex, hist_id, hist2_id, index, &
602	& evapot, evapot_corr, fluxsubli, lwabs, lwnet, netrad, transpir, &
603	& vevapflo, vevapnu, vevapp, vevapsno, vevapwet)
604	ENDIF
605
606	END SUBROUTINE hydraulic_arch_main
607
608
609	!! ================================================================================================================================
610	!! SUBROUTINE : hydraulic_arch_calc
611	!!
612	!>\BRIEF Calculate the amount of water that is available for the leaves for transpiration,
613	! accounting for effects of hydraulic architecture following (Hickler et al. 2006)
614	!!
615	!!\n DESCRIPTION : (1) Calculate soil water potential from soil water content as calculated in hydrol.f90
616	!! according to van Genuchten (1980). Soil water potential is weighted by the relative share of the
617	!! roots in each layer.
618	!! (2) Calculate pressure difference between leaves and soil (Whitehead 1998)
619	!! (3) Calculate root (Weatherly 1998), sapwood (Magnani et al. 2000, Tyree et al. 1994, Sperry et al. 1998)
620	!! and leave resistances (Sitch et al. 2003).
621	!! (4) Calculate the supply of water for transpiration by the leaves using Darcy's law (Whitehead 1998)
622	!!
623	!! RECENT CHANGE(S): None
624	!!
625	!! MAIN OUTPUT VARIABLE(S): :: transpir_supply
626	!!
627	!! REFERENCE(S) : +++COMPLETE+++
628	!! Hickler et al. 2006
629	!! van Genuchten (1980)
630	!! Whitehead 1998
631	!! Weatherly 1998
632	!! Magnani et al. 2000
633	!! Tyree et al. 1994
634	!! Sperry et al. 1998
635	!! Sitch et al. 2003
636	!!
637	!! FLOWCHART :
638	!!
639	!_ ================================================================================================================================
640
641	SUBROUTINE hydraulic_arch_calc (kjpindex, circ_class_biomass, circ_class_n, transpir, stempdiag, &
642	temp_air, swc, transpir_supply, njsc, veget, &
643	transpir_supply_column, ksoil, e_frac, vessel_loss, root_profile)
644
645
646	!! 0. Variable declaration
647
648	!! 0.1 Input variables
649	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
650	INTEGER(i_std), DIMENSION(:), INTENT (in) :: njsc !! Index of the dominant soil textural class
651	!! in the grid cell used for hydrology
652	!! (1-nscm, unitless)
653	REAL(r_std),DIMENSION(:,:),INTENT(in) :: stempdiag !! Soil temperature (K)
654	REAL(r_std),DIMENSION (:), INTENT (in) :: temp_air !! Air temperature (K)
655	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: swc !! Volumetric soil water content
656	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: circ_class_n !! number of trees within a circumference class (tree m-2)
657	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass components of the model tree
658	!! within a circumference class
659	!! class @tex $(g C ind^{-1})$ @endtex
660	REAL(r_std),DIMENSION (:,:), INTENT (in) :: veget !! Fraction of the vegetation that covers the soil
661	!! also project canopy cover - calculated by Pgap
662	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: ksoil !! Hydraulic conductivity mm/d
663	REAL(r_std), DIMENSION(:,:), INTENT (in) :: transpir !! Unstressed transpiration mm/dt
664	REAL(r_std),DIMENSION (:,:,:), INTENT(in) :: root_profile !! Normalized root mass/length fraction in each soil layer
665	!! (0-1, unitless)
666
667	!! 0.2 Output variables
668	REAL(r_std), DIMENSION(:,:), INTENT (out) :: transpir_supply !! Supply of water for transpiration. Same unit as transpir.
669	!! @tex $(mm m^{-2} leaf dt^{-1})$ @endtex
670	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: e_frac !! Fraction of water transpired supplied by individual layers (no units),
671	!! dimension(kjpindex,nvm,nslm,nstm)
672	REAL(r_std), DIMENSION(:,:), INTENT(out) :: vessel_loss !! Conductivity lost due to cavitation in the xylem (no unit).
673
674	!! 0.3 Modified variables
675
676
677	!! 0.4 Local variables
678	INTEGER(i_std) :: ier !! Error handling
679	INTEGER(i_std) :: ipts,ivm,icirc !! indices
680	INTEGER(i_std) :: ibdl,istm,jst !! indices
681	INTEGER(i_std) :: tmp_ncirc !! indices
682	REAL(r_std) :: R_root !! Hydraulic resistance of the roots (MPa s m-3)
683	REAL(r_std) :: R_leaf !! Hydraulic resistance of the leaves (MPa s m-3)
684	REAL(r_std) :: R_sap !! Hydraulic resistance of sapwood (MPa s m-3)
685	REAL(r_std) :: R_shoot !! Hydraulic resistance of shoot corrected for temperature (??)
686	REAL(r_std), DIMENSION(ncirc) :: height !! Vegetation height per circ_class_biomass (m)
687	REAL(r_std), DIMENSION(kjpindex,nslm,nstm) :: psi_soil !! Soil water potential per soil layer per soiltile higher
688	!! than -5 MPa (threshold for hygroscopic water) (MPa)
689	REAL(r_std), DIMENSION(ncirc) :: circ_class_transpir_supply !! Supply of water for transpiration per model tree
690	!! in each circumference class @tex $(m s^{-1}tree^{-1})$ @endtex
691	REAL(r_std) :: delta_P !! Pressure difference between leaves and soil
692	REAL(r_std), DIMENSION(nslm+1) :: z_soil !! Variable to store depth of the different
693	!! soil layers (m)
694	REAL(r_std), DIMENSION(nvm) :: k_sap_cav !! sapwood conductivity decreased by cavitation (m2 s-1MPa)
695	REAL(r_std), DIMENSION(nvm) :: rpc !! scaling factor for integral (unitless)
696	REAL(r_std), DIMENSION(kjpindex,nslm,nstm) :: psi_soiltile !! soil water potentential per soil layer per soil tile (MPa)
697	REAL(r_std), DIMENSION (kjpindex,nslm,nstm) :: swc_adj !! Volumetric soil water content, adjusted when it is lower
698	!! than or equal to the residual soil water content
699	REAL(r_std), DIMENSION(kjpindex,nvm) :: psi_soilroot !! Sum of soil water potential in each soil layer weighted
700	!! by the share of the roots in each layer (MPa)
701	INTEGER(i_std) :: i !! counter
702	REAL(r_std), DIMENSION(nlevels_tot, ncirc) :: cent_box !! Height of the mid-point of each box in the column that is made up
703	!! of the Pgap points
704	REAL(r_std), DIMENSION(nlevels_tot) :: delta_P_column !! Pressure difference between leaves and soil
705	REAL(r_std), DIMENSION(nlevels_tot,ncirc) :: circ_class_transpir_supply_column !! Supply of water for transpiration per model tree
706	!! in each circumference class @tex $(m s^{-1}tree^{-1})$ @endtex
707
708	REAL(r_std), DIMENSION(nlevels_tot,kjpindex,nvm):: transpir_supply_column !! Supply of water for transpiration
709	!! @tex $(mm m{-2} ground dt^{-1})$ @endtex
710	REAL(r_std), DIMENSION(kjpindex,nslm,nstm) :: e_max !! Maximum amount of available water for transpiration
711	!! per layer (mmol m2 s-1)
712	REAL(r_std), DIMENSION(kjpindex,nslm,nstm) :: res_soil !! Soil to root resistance per layer (s m-2)
713	REAL(r_std), DIMENSION(kjpindex,nslm,nstm) :: ksoilUNIT !! Hydraulic conductivity (from hydrol.90)
714
715	REAL(r_std), DIMENSION(nslm) :: lr !! Root length (m)
716	REAL(r_std), DIMENSION(nslm) :: rs !! Half the mean distance between roots (m)
717
718	REAL(r_std),DIMENSION(kjpindex,nvm,ncirc) :: R_root_circ !! Hydraulic resistance of the roots corrected for temperature
719	!! per circumference class (MPa s m-3))
720	REAL(r_std),DIMENSION(kjpindex,nvm,ncirc) :: R_shoot_circ !! Hydraulic resistance of the shoot corrected for temperature
721	!! per circumference class (MPa s m-3)
722	REAL(r_std),DIMENSION(kjpindex,nvm) :: R_root_tot !! Hydraulic resistance of the roots per pft
723	REAL(r_std),DIMENSION(kjpindex,nvm) :: R_shoot_tot !! Hydraulic resistance of the shoot per pft
724	REAL(r_std), DIMENSION(nslm) :: nroot_tmp !! Temporary variable to calculate the nroot
725	REAL(r_std), DIMENSION(kjpindex,nslm,nvm) :: tmp !! Temporary variable to write the root profile to xios
726
727	!_ =================================================================================================================================
728
729	!! Initialize local printlev
730	IF (l_first_sechiba) THEN
731	printlev_loc=get_printlev('hydrolarch')
732	ENDIF
733
734	! initialise some variables for the transpir_supply_column loop
735	i = 0
736	cent_box = zero
737	delta_P_column = zero
738	circ_class_transpir_supply_column = zero
739	transpir_supply_column = zero
740
741	! Save the hydraulic resitances
742	R_root_circ(:,:,:) = zero
743	R_shoot_circ(:,:,:) = zero
744	R_root_tot(:,:) = zero
745	R_shoot_tot(:,:) = zero
746
747	! To avoid that the plants will not experience water stress in the first
748	! day of the simulation. Transpir_supply is calculated before
749	! there is any biomass available.
750	transpir_supply(:,:) = val_exp
751	psi_soiltile(:,:,:) = -un
752	e_frac(:,:,:,:) = zero
753
754	! Initialize vessel_loss so it always has a value.
755	IF (hack_vessel_loss.LT.zero .OR. &
756	hack_vessel_loss.GT.un) THEN
757	! Set to zero. This is what we want for real simulations
758	vessel_loss(:,:) = zero
759	ELSE
760	! Use a constant fix value instead. This option
761	! was added for debugging and checking the sanity
762	! of the code
763	vessel_loss(:,:) = hack_vessel_loss
764	ENDIF
765
766	! Distinguish between the two soil map possibilities for these variables
767	! First we need to make sure there are allocated.
768
769	IF (.NOT. ALLOCATED(nvan)) THEN
770	ALLOCATE (nvan(nscm),stat=ier)
771	IF (ier /= 0) CALL ipslerr_p(3,'hydraulic_arch_calc','Pb in allocate of variable nvan','','')
772	ENDIF
773
774	IF (.NOT. ALLOCATED(avan)) THEN
775	ALLOCATE (avan(nscm),stat=ier)
776	IF (ier /= 0) CALL ipslerr_p(3,'hydraulic_arch_calc','Pb in allocate of variable avan','','')
777	ENDIF
778
779	IF (.NOT. ALLOCATED(mcr)) THEN
780	ALLOCATE (mcr(nscm),stat=ier)
781	IF (ier /= 0) CALL ipslerr_p(3,'hydraulic_arch_calc','Pb in allocate of variable mcr','','')
782	ENDIF
783
784	IF (.NOT. ALLOCATED(mcs)) THEN
785	ALLOCATE (mcs(nscm),stat=ier)
786	IF (ier /= 0) CALL ipslerr_p(3,'hydraulic_arch_calc','Pb in allocate of variable mcs','','')
787	ENDIF
788
789
790	nvan(:) = nvan_usda(:)
791	avan(:) = avan_usda(:)
792	mcr(:) = mcr_usda(:)
793	mcs(:) = mcs_usda(:)
794
795	!! 1. Calculate plant soil water potential
796
797	! Soil water potential for each soiltile in each layer from
798	! soil water content according to Van Genuchten (1980)
799
800	DO ipts = 1,kjpindex ! (loopref: ha1)
801
802	! Debug
803	IF (printlev_loc>=4) THEN
804	WRITE(numout,*) 'soil texture',njsc
805	WRITE(numout,*) 'swc', swc(ipts,:,:)
806	ENDIF
807	!-
808
809	!! 2.1 Calculate soil water potential
810	! Swc calculated in hydrol.f90 cannot be lower than or equal
811	! to the residual soil water content for each soiltile.
812	DO ibdl = 1, nslm ! (loopref: ha2)
813
814	DO istm = 1,nstm ! (loopref: ha3)
815
816	IF (swc(ipts,ibdl,istm) .LE. mcr(njsc(ipts))) THEN
817
818	!+++CHECK+++
819	!Explain why you do this
820	swc_adj(ipts,ibdl,istm) = mcr(njsc(ipts)) + .00001
821	!+++++++++++
822
823	ELSE
824
825	swc_adj(ipts,ibdl,istm) = swc(ipts,ibdl,istm)
826
827	ENDIF
828
829	ENDDO ! istm = 1,nstm (loopref: ha3)
830
831	! Soil water potential for each soiltile in each layer from
832	! soil water content according to Van Genuchten (1980). Convert
833	! unit avan from mm-1 to m-1. Convert units soil matric potential
834	! from m to MPa by multiplying with cte_grav and rho_h2o.
835	! Soiltile = fraction of each soil column
836	! (three soil columns=bare soil, tree, grass)
837	DO istm = 1,nstm ! (loopref: ha4)
838
839	IF (swc_adj(ipts,ibdl,istm) .GT. min_stomate .AND. &
840	(swc_adj(ipts,ibdl,istm) - mcr(njsc(ipts))) .GT. min_stomate) THEN
841
842
843	psi_soiltile(ipts,ibdl,istm) = &
844	- ( cte_grav * rho_h2o * kilo_to_unit * &
845	( 1/(avan(njsc(ipts))kilo_to_unit) ) &
846	( ( (swc_adj(ipts,ibdl,istm) - mcr(njsc(ipts)) ) / &
847	( mcs(njsc(ipts)) - mcr(njsc(ipts)) ) ) ** &
848	( -1 / (1-1/nvan(njsc(ipts))) ) - 1) ** &
849	(1/nvan(njsc(ipts)) ) ) / mega_to_unit
850
851	ELSE
852
853	psi_soiltile(ipts,ibdl,istm) = zero
854
855	ENDIF
856
857	ENDDO ! DO istm = 1,nstm (loopref: ha4)
858
859	! Threshold for soil water potential: it cannot be lower than the
860	! soil water potential for hygroscopic water (= -5 MPa, Larcher 2001)
861	psi_soil(ipts,ibdl,:) = MAX(psi_soiltile(ipts,ibdl,:),-5.)
862
863
864	ENDDO ! DO ibdl = 1, nslm ! (loopref: ha2)
865
866	! following initialisation of swc
867	l_first_sechiba = .FALSE.
868
869	! Debug
870	IF (printlev_loc>=4) THEN
871	WRITE(numout,*) 'psi_soiltile', psi_soiltile(ipts,:,:)
872	WRITE(numout,*) 'psi_soil', psi_soil
873	ENDIF
874	!-
875
876	!! 3. Calculate soil to root and plant resistances
877
878	!! 3.1 Calculate root profile
879
880	! This is now being done in hydrol.f90. Given hydraulic_arch is called
881	! before hydrol.f90, the root_profile runs one time step behind the
882	! soil moisture profile. This does not appear as an important issue.
883	! In case the static root profile is used the root profile is fixed and
884	! there is no issue at all. In case a dynamic root profile is used,
885	! this 30 minute delay would occur. See hydrol_dynamic_root_profile
886	! for some details/thinking on the issue.
887
888	!! 3.2 Calculate water transport from the soil-root zone up to the leaves
889	! Initialize psi_soilroot to avoid problems with sending unitialized values
890	! to xios
891	psi_soilroot(ipts,1) = xios_default_val
892
893	DO ivm = 2,nvm ! (loopref: ha6)
894
895	!! Calculate PFT-specific characteristics
896	IF (is_tree(ivm)) THEN
897
898	! Calculate for each circumference classes
899	tmp_ncirc = ncirc
900
901	! Calculate the effective height (m) from biomass. The effective
902	! value is used here because the water needs to be transported
903	! through both the belowground and aboveground biomass. Using the
904	! height (thus not the effective height) would only account for the
905	! height due to the aboveground biomass.
906	height(:) = wood_to_height_eff(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
907
908	ELSE ! is_tree(ivm)
909
910	! There are no circumference classes for grasses
911	! and crops. Biomasses are stored in first circumference
912	! class
913	tmp_ncirc = 1
914
915	! Calculate height (m) from biomass. Use the function
916	! biomass_to_lai to account for dynamic sla
917	height(1) = biomass_to_lai(&
918	circ_class_biomass(ipts,ivm,1,ileaf,icarbon),ivm) * lai_to_height(ivm)
919	ENDIF
920
921	!! 3.2.1 Calculate soil to root resistance
922	! Calculate soil water potential weighted for root density
923	! Note that the second soiltile in swc is for trees.
924	! This implementation relies on a strong assumptions. First,
925	! roots are assumed to follow a rigid distribution irrespective
926	! of how the water is distributed in the soil. So even if on
927	! average lots of water is available in deeper layer, the modeled
928	! plants will grow most roots in the top layers. A detailed description
929	! allows for an adjustment of the soil water potential (psi_soilroot)
930	! due to soil to root resistance. Psi_soil is weighted by the fraction of
931	! evapotranspiration supplied by each layer (E_frac), which depends on
932	! the soil to root resistance per layer (res_soil) accounting for
933	! different root properties and hydraulic conductivity.
934
935	!Select soil tile for this pft
936	jst = pref_soil_veg(ivm)
937
938	! We can only make these calculations if there is actual root
939	! biomass.
940	IF(SUM(circ_class_biomass(ipts,ivm,:,iroot,icarbon)) .GE. min_stomate ) THEN
941
942	DO ibdl = 1, nslm
943
944	! Calculate root length and one-half the mean distance between roots (rs)
945	! following Bonan et al. (2014) GMD. eq A23
946
947	! Calculate root length (m) using parameter specific root length (srl).
948	! Convert gC/ind to g/ind (deux).
949	lr(ibdl)=root_profile(ipts,ivm,ibdl)deux &
950	SUM(circ_class_biomass(ipts,ivm,:,iroot,icarbon)*&
951	circ_class_n(ipts,ivm,:)) * srl(ivm)
952
953	! Grasses and crops have only roots to one meter deep
954	! hence lr(ibdl) will be zero fo the deeper soil layers
955	IF (lr(ibdl).GT.min_stomate) THEN
956	rs(ibdl) = (un/(pilr(ibdl)))*(undemi)
957	ELSE
958	rs(ibdl) = zero
959	END IF
960
961	! Convert units for the hydraulic conductivity from mm/d to m/s
962	ksoilUNIT(ipts,ibdl,jst)=ksoil(ipts,ibdl,jst)/(mm_m * one_day)
963
964	! The soil to root resistance in each layer (res_soil) is calculated before
965	! the contribution from each layer to the evapotransporation e_frac.
966	! froot gives the mean radius of fine roots (m).
967	! Unit of res_soil is s m-2. Further
968	! unit conversion is not needed because e_max is normalized.
969	IF (rs(ibdl).GT.min_stomate) THEN
970	res_soil(ipts,ibdl,jst)=(LOG(rs(ibdl)/r_froot(ivm)))/&
971	(deuxpilr(ibdl)*ksoilUNIT(ipts,ibdl,jst))
972	ELSE
973	! If there are no roots in a specific layer the resistance is
974	! infinite.
975	res_soil(ipts,ibdl,jst)=large_value
976	END IF
977
978	! Psi_root is a poorly measured parameter. Since psi_leaf_min is
979	! frequently measured and a psi_root that is lower than psi_leaf
980	! will not result in more more water transport to the leaves anyway, we
981	! prefer to set psi_root equal to psi_leaf rather than choosing an
982	! arbitrary value that is unrelated to to the prescribed psi_leaf.
983	psi_root(ivm)=psi_leaf(ivm)
984	e_max(ipts,ibdl,jst)= (psi_soil(ipts,ibdl,jst)-psi_root(ivm))/&
985	res_soil(ipts,ibdl,jst)
986
987	! If we had to divide by an infinite resistance (where infinite was
988	! represented by large_value (=10e+33), e_max should less than
989	! min_sechiba and will be set to zero.
990	IF (e_max(ipts,ibdl,jst) .LT. min_sechiba) e_max(ipts,ibdl,jst) = zero
991
992	END DO
993
994	! If there is a contribution to evapotransporation, the fraction
995	! coming from each layer is calculated first, followed by calculations
996	! of the weighted soilroot water potential (psi_soilroot). When there is no
997	! contribution, psi_soilroot is set very low.
998	IF( SUM(e_max(ipts,:,jst)) .GT. min_sechiba) THEN
999
1000	!+++CHECK+++
1001	! There is no need to have the nst dimension in e_frac because each PFT
1002	! is linked to just a single soil tile. Could be simplified.
1003	DO ibdl = 1, nslm
1004	e_frac(ipts,ivm,ibdl,jst)=e_max(ipts,ibdl,jst)/SUM(e_max(ipts,:,jst))
1005	ENDDO
1006	!+++++++++++
1007
1008	IF (hack_e_frac) THEN
1009	! Simplify the calculation of psi_soilroot by not accounting
1010	! for root length. Overwrite e_frac with root_profile.
1011	e_frac(ipts,ivm,:,jst) = root_profile(ipts,ivm,:)
1012	END IF
1013
1014	! Calculate psi_soilroot
1015	psi_soilroot(ipts,ivm) = SUM( psi_soil(ipts,:,jst)*e_frac(ipts,ivm,:,jst))
1016
1017	ELSE
1018
1019	! The calculation of psi_soil is not realistic because psi_root > psi_soil
1020	! Use a suboptimal solution by prescribing e_frac and psi_soilroot
1021	e_frac(ipts,ivm,:,jst)= un/nslm
1022	psi_soilroot(ipts,ivm)=-4
1023
1024	END IF
1025
1026	! Debug
1027	IF(printlev_loc .GE. 4 .AND. ivm .eq. test_pft) THEN
1028	WRITE(numout,*)'hydraulic_arch_calc in biomass statement, ivm, jst',ivm, jst
1029	WRITE(numout,*)'psi_soil', psi_soil(ipts,:,jst)
1030	WRITE(numout,*)'lr',lr(:)
1031	WRITE(numout,*)'rs', rs(:)
1032	WRITE(numout,*)'ksoilUNIT', ksoilUNIT(ipts,:,jst)
1033	WRITE(numout,*)'res_soil',res_soil(ipts,:,jst)
1034	WRITE(numout,*)'e_max',e_max(ipts,:,jst)
1035	WRITE(numout,*)'e_frac', e_frac(ipts,ivm,:,jst)
1036	WRITE(numout,*)'psi_soilroot',psi_soilroot(ipts,ivm)
1037	ENDIF
1038
1039	ELSE ! IF(biomass .gt. min_stomate) THEN
1040
1041	! Only do the calculations for biomass
1042	e_frac(ipts,ivm,:,jst)= zero
1043	psi_soilroot(ipts,ivm) = SUM(psi_soil(ipts,:,jst)*e_frac(ipts,ivm,:,jst))
1044
1045	ENDIF ! (biomass .GT. min_stomate) THEN
1046
1047
1048	!! 3.2.2 Calculate plant resistances
1049	! Calculate plant resistance to water transport according to the principle
1050	! of an electric circuit. The code is largely based on Hickler et al. 2006.
1051
1052	! Set transpir_supply to zero whem transpir is zero (e.g. at night),
1053	! so their daily means are comparable.
1054	IF ( SUM(circ_class_biomass(ipts,ivm,:,ileaf,icarbon)) .GT. min_sechiba .AND. &
1055	SUM(circ_class_biomass(ipts,ivm,:,iroot,icarbon)) .GT. min_sechiba .AND. &
1056	transpir(ipts,ivm) .GT. min_sechiba ) THEN
1057
1058	DO icirc = 1,tmp_ncirc ! (loopref ha7)
1059
1060	! This ensures that the loop goes over all circ_classes for trees and only
1061	! the first circ_class for grasses and crops
1062	IF (circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon) .GT. min_sechiba .AND. &
1063	circ_class_biomass(ipts,ivm,icirc,iroot,icarbon) .GT. min_sechiba ) THEN
1064
1065	!! 3.2.2.1 Root resistance
1066	! Root resistance is in the first order related to root mass
1067	! and expression was proposed by Weatherly, 1982). Convert
1068	! units of circ_class_biomss i.e. gC/ind to kg/ind
1069	! Unit of R_root is (MPa s m-3 ind-1)
1070	R_root = un / (k_root(ivm) * 2 * &
1071	circ_class_biomass(ipts,ivm,icirc,iroot,icarbon) &
1072	/ kilo_to_unit)
1073
1074	!! 3.2.2.2 Sapwood resistance
1075	! Calculate k_sap_cav as a function of psi_soilroot and the pft specific parameter
1076	! psi_50 which representing its vulnerability to cavitation, following
1077	! Tyree et al. (1994) and Sperry et al. (1998).
1078	IF (is_tree(ivm)) THEN
1079
1080	!! !+++++CHECK: psi_50 should be lower than 0 to have cavitation, so it seems we never
1081	!! reduce k_sap due to cavitation for the moment!
1082	IF( (psi_soilroot(ipts,ivm) == zero) .AND. (psi_50(ivm) .LE. zero) )THEN
1083
1084	! This seems to cause a strange numerical issue, even though we should
1085	! just be taking the exponential of zero (i.e., one).
1086	k_sap_cav(ivm) = k_sap(ivm)
1087
1088	ELSE
1089	!++++ TEMP +++++
1090	IF((psi_soilroot(ipts,ivm)/psi_50(ivm)) .LT. zero)THEN
1091	!++++++ CHECK +++++++++
1092	! If this is true, we generally try to take this ratio to
1093	! an odd power, which results in an imaginary number. So
1094	! we use this simple solution for the moment.
1095	! WRITE(numout,*) 'psi_soilroot', psi_soilroot(ipts,ivm)
1096	! WRITE(numout,*) 'ipts', ipts
1097	! WRITE(numout,*) 'ivm', ivm
1098	! CALL ipslerr_p (3,'hydraulic_arch_calc', &
1099	! 'this quantity should never be negative', &
1100	! 'See the output file for more details.','')
1101	! A temporary fix. We cannot have the stop statement anymore,
1102	! since capping psi_soilroot at 0.0 above results in many PFTs
1103	! dying. This needs to be addressed.
1104	k_sap_cav(ivm) = k_sap(ivm)
1105	ELSE
1106	k_sap_cav(ivm) = k_sap(ivm) * &
1107	EXP(-((psi_soilroot(ipts,ivm)/psi_50(ivm)) ** &
1108	c_cavitation(ivm)))
1109	ENDIF
1110
1111	ENDIF
1112
1113	! If k_sap_cav is too low, we can get a divide by zero
1114	! error in the next couple lines or vessel_loss could exceed
1115	! one. So this is a protection against that.
1116	IF(k_sap_cav(ivm) .LE. min_sechiba)THEN
1117
1118	k_sap_cav(ivm) = min_sechiba
1119
1120	ENDIF
1121
1122	! Caculation of the PLC (percentage of lost conductivity
1123	! in the xylem).
1124	IF (ok_vessel_mortality .AND. &
1125	(hack_vessel_loss.LT.zero .OR. &
1126	hack_vessel_loss.GT.un)) THEN
1127
1128	! The ratio between k_sap_cav and k_sap
1129	! gives the fraction of functional vessels remaining. The
1130	! PLC will be used in module stomate_turnover.f90 to
1131	! compute mortality.
1132	vessel_loss(ipts,ivm) = un - (k_sap_cav(ivm) / k_sap(ivm))
1133
1134	END IF
1135
1136	! Following the expression by Magnani et al 2000, sapwood resistance
1137	! is given as a function of sapwood area. Convert Sapwood mass to
1138	! sapwood area.
1139	! Units of R_sap is (MPa s m-3 ind-1).
1140	R_sap = (height(icirc)*2 tree_ff(ivm)pipe_density(ivm)) / (k_sap_cav(ivm) &
1141	( circ_class_biomass(ipts,ivm,icirc,isapabove,icarbon) + &
1142	circ_class_biomass(ipts,ivm,icirc,isapbelow,icarbon)))
1143
1144	ELSE ! IF (is_tree(ivm)) THEN
1145
1146	!+++CHECK+++
1147	! We have about the same leaf and root restistances and masses
1148	! as for trees but lack a resistance in the sapwood. This is resulting
1149	! in a huge supply because we lack part of the resistance. This was
1150	! the motivation to create an R_sap for grasses and roots. This equation
1151	! is made up and should be carefully checked against physiological
1152	! studies for grasses and crops on the issue.
1153	R_sap = (height(icirc)*2 pipe_density(ivm)) / (k_sap(ivm) * &
1154	circ_class_biomass(ipts,ivm,icirc,isapabove,icarbon))
1155	!+++++++++++
1156
1157	ENDIF ! IF (is_tree(ivm)) THEN
1158
1159	!! 3.2.2.3 Leaf resistance
1160	! Leaf resistance related to foliar projective area (Hickler 2006) or
1161	! LAI (more consistent with pipe model)? First fpc= 1-e*(-0.5LAI) was used
1162	! instead of lai because it is assumed that only sunlit leaves contribute to transpiration.
1163	! Then it was replaced by veget which in itself is calculated in Pgap and thus accounts for
1164	! LAI and canopy structure. This was still inconsistent with the calculation of the
1165	! atmospheric demand for transpiration. To improve consistency the coupled
1166	! lai is now used as done in the energy budget calculations.
1167	! R_leaf = 1/(k_leaf * fpc)
1168	!!$ R_leaf = un / (k_leaf(ivm) * &
1169	!!$ (un - EXP(-0.5 * circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon) * &
1170	!!$ sla(ivm))))
1171	! Units of R_leaf is (MPa s m-3), because biomass_to_couple_lai
1172	! takes the biomass of the leaves at the tree level, and
1173	! hence we get the leaf area per individual, which are in m2.
1174	R_leaf = un / (k_leaf(ivm) * &
1175	biomass_to_coupled_lai(circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon), &
1176	veget(ipts,ivm),ivm))
1177
1178	!! 3.2.2.4 Temperature dependance of resistance
1179	! Decrease root resistance with soil temperature as a result of
1180	! changes in water viscosity (Cochard et al 2000). grnd_80 is the
1181	! index of the soil layer that is 80 cm deep. This 80 cm is arbitrary
1182	! and was used because we already used the same variable in wind throw
1183	! so the variable was already calculated.
1184	! Decrease shoot resistance with air temperature (Cochard et al. 2000)
1185	R_root_circ(ipts,ivm,icirc) = R_root / (a_viscosity(1) + a_viscosity(2) * &
1186	(stempdiag(ipts,grnd_80) - zeroCelsius))
1187	R_shoot_circ(ipts,ivm,icirc) = (R_leaf + R_sap) / (a_viscosity(1) + a_viscosity(2) * &
1188	(temp_air(ipts) - zeroCelsius))
1189
1190	! (JR 150115) in this version we try to calculate delta_P against height
1191	! calculate the 'mid-point' of the P_gap points
1192
1193	IF (.NOT. ok_impose_can_structure ) THEN
1194	!DO i = 1, nlevels_tot - 1
1195	! cent_box(i, icirc) = z_array_out3(1,ivm,1,i) + (z_array_out3(1,ivm,1,i+1) - z_array_out3(1,ivm,1,i)) / 2.
1196	!END DO
1197	!cent_box(nlevels_tot, icirc) = z_array_out3(1,ivm,1,nlevels_tot) + &
1198	! & ( h_array_out3(1,ivm,1, nlevels_tot) - z_array_out3(1,ivm,1,nlevels_tot) ) / 2.0
1199	! DO i = 1, nlevels_tot
1200	! WRITE(numout, *) 'cent_box(i, icirc) is: ', cent_box(i, icirc)
1201	! END DO
1202	! DO i = 1, nlevels_tot
1203	! delta_P_column(i) = psi_leaf(ivm) - psi_soilroot(ipts,ivm) + &
1204	! (cent_box(i, icirc)cte_gravrho_h2o) / kilo_to_unit
1205	! END DO
1206	END IF ! (ok_impose_can_structure ) THEN
1207
1208	!! 3.2.3 Calculate pressure difference between leaves and soil
1209	! Defined by Whitehead 1998. Convert units to MPa
1210	! Transpir_supply is to a large extend determined by delta_P. DeltaP
1211	! in turn is sensitive to two prescribed parameters: psi_leaf_min
1212	! and psi_root. Psi_leaf_min is rather well documented but using it
1213	! implies that we always calculate the maximum transpir_supply which
1214	! results in drying the soil too quickly in dry regions
1215	! (CHECK:still true in cleaned version?). Psi_root can
1216	! be tuned (currently set to psi_leaf, see later) because it is poorly observed.
1217	! Some of the expected dynamics might not be captured, for example, when
1218	! a drought is emerging, transpiration increases and then suddenly collapses
1219	! (CHECK if this behaviour is still occurs).
1220	! In reality it is expected that the soil would dry and that a drier soil would reduce
1221	! the water the plant can transport. Owing to the static character of
1222	! this scheme, it looks like its limits are reached. Further improving
1223	! water stress would require a more dynamic scheme that calculates
1224	! psi_leaf_min and psi_root as a function of the plant and soil water status and respectively.
1225	delta_P = psi_leaf(ivm) - psi_soilroot(ipts,ivm) + &
1226	(height(icirc)cte_gravrho_h2o) / kilo_to_unit
1227
1228
1229	!! 3.2.4 Calculate water supply for transpiration.
1230	! First for a model tree in each circ_class (m/s/tree), then
1231	! for each pft.
1232
1233	! Using Dracy's law (Slatyer 1967, Whitehead 1998)
1234	circ_class_transpir_supply(icirc) = (-delta_P) / &
1235	(R_root_circ(ipts,ivm,icirc) + R_shoot_circ(ipts,ivm,icirc))
1236
1237
1238	IF (.NOT. ok_impose_can_structure ) THEN
1239	DO i = 1, nlevels_tot
1240	circ_class_transpir_supply_column(i, icirc) = (-delta_P_column(i)) / &
1241	(R_root_circ(ipts,ivm,icirc) + R_shoot_circ(ipts,ivm,icirc))
1242	END DO
1243	END IF ! (ok_impose_can_structure ) THEN
1244
1245	ELSE
1246
1247	!set transpir_supply to zero if there is no leaf or root biomass in the circ_class
1248	circ_class_transpir_supply(icirc) = zero
1249
1250	ENDIF
1251
1252	!---TEMP---
1253	IF (printlev_loc>=4 .AND. ivm .EQ. test_pft) THEN
1254	WRITE(numout,*) 'root_profile', root_profile(ipts,ivm,:)
1255	WRITE(numout,*) 'k_sap_cavitation, biomass(iroot)', k_sap_cav(ivm), &
1256	2 * circ_class_biomass(ipts,ivm,icirc,iroot,icarbon) / 1000.
1257	WRITE(numout,*) 'sapwood area, ', ( circ_class_biomass(ipts,ivm,icirc,isapabove,icarbon) + &
1258	circ_class_biomass(ipts,ivm,icirc,isapbelow,icarbon)) / (tree_ff(ivm) * pipe_density(ivm)) / &
1259	height(icirc)
1260	WRITE(numout,) 'fpc, ',(1-EXP(-0.5circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon)*sla(ivm)))
1261	WRITE(numout,*) 'R_root, R_sap, R_leaf', R_root_circ(ipts,ivm,icirc), R_sap, R_leaf
1262	WRITE(numout,*) 'R_shoot', R_shoot_circ(ipts,ivm,icirc)
1263	WRITE(numout,*) 'delta_P', delta_P
1264	ENDIF
1265
1266	ENDDO ! DO icirc = 1,tmp_ncirc (loopref ha7)
1267
1268
1269	! Calculating transpiration supply at different heights here:
1270	IF (.NOT. ok_impose_can_structure ) THEN
1271
1272	DO i = 1, nlevels_tot
1273	transpir_supply_column(i, ipts, ivm) = &
1274	(SUM(circ_class_transpir_supply_column(i, :) * circ_class_n(ipts,ivm,:))) * &
1275	dt_sechiba * kilo_to_unit
1276
1277	END DO ! i = 1, nlevels_tot
1278
1279
1280
1281	! Set transpir_supply to zero when negative.
1282	! (When psisoilroot becomes very low, delta_P becomes positive and transpir_supply negative)
1283	DO i = 1, nlevels_tot
1284	IF (transpir_supply_column(i, ipts, ivm) .LT. zero) THEN
1285	WRITE(numout,*) 'transpir_supply before set to zero when negative', transpir_supply(ipts,ivm), ivm
1286	WRITE(numout,*) 'psisoilroot', psi_soilroot(ipts,ivm)
1287	transpir_supply_column(i, ipts, ivm) = MAX(transpir_supply_column(i, ipts, ivm), zero)
1288	ELSE
1289	!do nothing
1290	ENDIF
1291	END DO ! i = 1, nlevels_tot
1292
1293	END IF ! (ok_impose_can_structure ) THEN
1294
1295
1296	! Calculate transpir_supply per pft. Unit conversion, convert m/s/tree to mm/m^2/half_hour.
1297	! Multiply with veget to have the same unit as transpir in sechiba,
1298	! transpir is calculated by leaf area not by land area.
1299	IF (is_tree(ivm)) THEN
1300	! transpir_supply per PFT (sum of diameter classes) per timestep
1301	transpir_supply(ipts,ivm) = &
1302	(SUM(circ_class_transpir_supply(:) * circ_class_n(ipts,ivm,:))) * &
1303	dt_sechiba * kilo_to_unit * veget(ipts,ivm)
1304	ELSE
1305	transpir_supply(ipts,ivm) = &
1306	(circ_class_transpir_supply(1) * circ_class_n(ipts,ivm,1)) * &
1307	dt_sechiba * kilo_to_unit* veget(ipts,ivm)
1308	ENDIF
1309
1310	!+++++++++++
1311	! calculate the hydraulic resistances per PFT
1312	R_root_tot(ipts,ivm) = (SUM(R_root_circ(ipts,ivm,:)*circ_class_n(ipts,ivm,:)))
1313	r_shoot_tot(ipts,ivm) = (SUM(R_shoot_circ(ipts,ivm,:)*circ_class_n(ipts,ivm,:)))
1314
1315	! Set transpir_supply to zero when negative.
1316	! (When psisoilroot becomes very low, delta_P becomes positive and transpir_supply negative)
1317	IF (transpir_supply(ipts,ivm) .LT. zero) then
1318	IF(printlev_loc.GT.4 .and. ivm .eq. test_pft)THEN
1319	WRITE(numout,*) 'transpir_supply before set to zero when negative', transpir_supply(ipts,ivm), ivm
1320	WRITE(numout,*) 'ipts, ivm:', ipts,ivm
1321	WRITE(numout,*) 'transpir_supply: ', transpir_supply(ipts,ivm)
1322	WRITE(numout,*) 'psisoilroot: ', psi_soilroot(ipts,ivm)
1323	ENDIF
1324	transpir_supply(ipts,ivm) = MAX(transpir_supply(ipts,ivm), zero)
1325	ELSE
1326	!do nothing
1327	ENDIF
1328
1329
1330	! Hydraulic architecure is not defined i.e. no biomass
1331	ELSE
1332
1333	transpir_supply(ipts,ivm) = zero
1334	psi_soilroot(ipts,ivm) = zero
1335
1336	! DO i = 1, nlevels_tot
1337	! transpir_supply_column(i, ipts, :) = zero
1338	! END DO ! i = 1, nlevels_tot
1339
1340
1341	ENDIF ! biomass leaves and roots.gt.zero
1342
1343	ENDDO ! ivm = 2,nvm (loopref: ha6)
1344
1345	ENDDO ! kjpindexipts ! DO ipts = 1,kjpindex (loopref: ha1)
1346
1347	! Debug
1348	IF (printlev_loc>=4) THEN
1349	WRITE(numout,*) 'transpir_supply, ', transpir_supply(:,:)
1350	ENDIF
1351	!-
1352
1353	! Send output variables to xios
1354	CALL xios_orchidee_send_field("psi_soil",psi_soil)
1355	CALL xios_orchidee_send_field("psi_soilroot",psi_soilroot)
1356	CALL xios_orchidee_send_field("transpir_unstressed",transpir*one_day/dt_sechiba)
1357	CALL xios_orchidee_send_field("transpir_supply",transpir_supply*one_day/dt_sechiba)
1358	CALL xios_orchidee_send_field("R_root",R_root_tot)
1359	CALL xios_orchidee_send_field("R_shoot",R_shoot_tot)
1360	CALL xios_orchidee_send_field("VESSEL_LOSS",vessel_loss)
1361
1362	END SUBROUTINE hydraulic_arch_calc
1363
1364
1365	!! ==============================================================================================================================
1366	!! SUBROUTINE : mleb_diffuco_trans_co2_short
1367	!!
1368	!>\BRIEF This subroutine recalculate carbon assimilation and stomatal
1369	!! conductance following vbeta3 recalculation in sechiba.f90
1370	!!
1371	!! DESCRIPTION :\n
1372	!! *** General:\n
1373	!! cresist, sum_rveget_rstruct, gs, assimi, and ci are recalculated sequentially according to
1374	!! the original equations in diffuco_trans_co2. Because this routine places inside of the loop
1375	!! for pixel and pft, the subroutine itself doesn't contain loops for pixel and pft but uses
1376	!! ji and jv as inputs. Every equation should be identical with diffuco_trans_co2, so any
1377	!! updates and changes from diffuco_trans_co2 should be applied here, too.
1378
1379	!_ ================================================================================================================================
1380
1381	SUBROUTINE mleb_diffuco_trans_co2_short(assimi_lev, assim_param, Ca, gstot_frac, ji, jv, kjpindex, &
1382	Light_Abs_Tot, Light_Tran_Tot, lai, lai_per_level, pb, qair, q_cdrag, &
1383	qsintmax, qsintveg, ratio_vbeta3, rstruct, rveget, swdown, temp_air, temp_growth, u, &
1384	v, vbeta23, vbeta3, veget_max, veget, info_limitphoto, JJ_out, cimean, cresist, &
1385	new_leaf_ci, vbeta3pot, gsmean, gpp)
1386	!
1387	!! 0. Variable and parameter declaration
1388	!
1389	!
1390	!! 0.1 Input variables
1391	!
1392	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: assimi_lev
1393	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: assim_param !! min+max+opt temps (K), vcmax, vjmax
1394	!! for photosynthesis ($\mumolm^{-2}s^{-1}$)
1395	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: Ca !! CO2 concentration inside the canopy
1396	!! @tex ($\mu mol mol^{-1}$) @endtex
1397	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: gstot_frac !! fraction of gs over nlevel_tot (-)
1398	INTEGER(i_std), INTENT(in) :: ji !! index for pixel
1399	INTEGER(i_std), INTENT(in) :: jv !! index for PFT
1400	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1401	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: Light_Abs_Tot !! Absorbed radiation per level for photosynthesis
1402	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: Light_Tran_Tot !! Transmitted radiation per level for photosynthesis
1403	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Surface foliaire
1404	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lai_per_level !! This is the LAI per vertical level
1405	!! @tex $(m^{2} m^{-2})$
1406	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Lowest level pressure (hPa)
1407	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
1408	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: q_cdrag !! Surface drag coefficient (-)
1409	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: qsintmax !! Maximum water on vegetation
1410	!! @tex ($kg m^{-2}$) @endtex
1411	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: qsintveg !! Water on vegetation due to interception
1412	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: ratio_vbeta3 !! Reduction rate of vbeta3 when transpir_supply
1413	!! is smaller than demand
1414	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: rstruct !! structural resistance @tex ($s m^{-1}$) @endtex
1415	!! @tex ($\mu mol mol^{-1}$) @endtex
1416	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: rveget !! stomatal resistance of vegetation
1417	REAL(r_std), DIMENSION(:), INTENT(in) :: swdown !! Downwelling short wave flux
1418	!! @tex ($W m^{-2}$) @endtex
1419	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: temp_air !! Air temperature at first atmospheric model layer (K)
1420	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
1421	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: u !! Lowest level wind speed @tex ($m s^{-1}$) @endtex
1422	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: v !! Lowest level wind speed @tex ($m s^{-1}$) @endtex
1423	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: vbeta23 !! Beta for fraction of wetted foliage that will
1424	!! transpire (unitless)
1425	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: vbeta3 !! Beta for Transpiration (unitless)
1426	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! Maximum vegetation fraction of each PFT inside
1427	!! the grid box (0-1, unitless)
1428	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation type (-)
1429	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot+1), INTENT(in):: info_limitphoto!!Save information of limitation for photosynthesis
1430	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in) :: JJ_out
1431
1432	!! 0.2 Output variables
1433	!
1434	!! 0.3 Modified variables
1435	!
1436	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: cimean !! mean intercellular CO2 concentration
1437	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: cresist !! coefficient for resistances (-)
1438	REAL(r_std), DIMENSION(kjpindex,nvm,nlai), INTENT(inout) :: new_leaf_ci !! intercellular CO2 concentration (ppm)
1439	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: vbeta3pot !! Beta for Potential Transpiration
1440	!! @tex ($s m^{-1}$) @endtex
1441	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: gsmean !! mean stomatal conductance to CO2 (umol m-2 s-1)
1442	REAL(r_Std), DIMENSION(kjpindex,nvm), INTENT(inout) :: gpp !! Assimilation ((gC m^{-2} s^{-1}), total area)
1443
1444	!! 0.4 Local variables
1445	!
1446	REAL(r_std) :: assim_coeff_a, assim_coeff_b, &
1447	assim_coeff_x1, assim_coeff_x2, assim_coeff_x3 !! Coefficient based on analytic solution in diffuco
1448	REAL(r_std), DIMENSION(nlevels_tot) :: assimi !! assimilation (at a specific LAI level)
1449	REAL(r_std) :: assimtot !! total assimilation
1450	REAL(r_std) :: Cs_star !! Cs -based CO2 compensation point
1451	!! in the absence of Rd (ubar)
1452	REAL(r_std) :: ci_star !! Ci -based CO2 compensation point in the absence of Rd (ubar)
1453	REAL(r_std) :: fcyc !! Fraction of electrons at PSI that follow cyclic
1454	!! transport around PSI (-)
1455	REAL(r_std) :: fvpd !! Factor for describing the effect of leaf-to-air vapour
1456	!! difference on gs (-)
1457	REAL(r_std) :: gamma_star !! CO2 compensation point (ppm)
1458	REAL(r_std) :: gm !! Mesophyll diffusion conductance
1459	!! (molCO2 mÃ¢2 sÃ¢1 barÃ¢1)
1460	REAL(r_std) :: gb !! Boundary-layer conductance (mol m^{-2} s^{-1} bar{-1})
1461	REAL(r_std), DIMENSION(nlevels_tot) :: gs !! stomatal conductance to CO2
1462	!! @tex ($\mol m^{-2} s^{-1}$) @endtex
1463	REAL(r_std) :: gstot !! total stomatal conductance to H2O
1464	REAL(r_std) :: g0var
1465	INTEGER(i_std) :: ilev !! index
1466	REAL(r_std) :: J2 !! Rate of all e-transport through PSII
1467	!! (umol e- m-2 s-1)
1468	REAL(r_std) :: jwind !! Wind module (m s^{-1})
1469	REAL(r_std) :: Kmc !! MichaelisâMenten constant of Rubisco for CO2 (mubar)
1470	REAL(r_std) :: KmO !! MichaelisâMenten constant of Rubisco for O2 (mubar)
1471	REAL(r_std) :: lai_total !! total LAI for column in question
1472	REAL(r_std) :: leafN_top !! Leaf N content - top of canopy (gN m-2[leaf])
1473	REAL(r_std), DIMENSION(nlevels_tot) :: light_tran_to_level !! Cumulative amount of light transmitted
1474	!! to a given level
1475	REAL(r_std) :: low_gamma_star !! Half of the reciprocal of Sc/o (bar bar-1)
1476	REAL(r_std) :: N_profile !! N-content reduction factor for a specific canopy depth (unitless)
1477	REAL(r_std) :: N_Vcmax !! Nitrogen level dependance of Vcmacx and Jmax
1478	REAL(r_std) :: ntot !! Total canopy (eg leaf) N (gN m-2[ground])
1479	REAL(r_std) :: nue !! Nitrogen use Efficiency with impact of leaf
1480	REAL(r_std) :: Obs !! Bundle-sheath oxygen partial pressure (ubar)
1481	REAL(r_std), DIMENSION(kjpindex) :: qsatt !! surface saturated humidity at 2m (??)
1482	REAL(r_std), DIMENSION(nlevels_tot) :: Rd !! Day respiration (respiratory CO2 release other than
1483	!! by photorespiration) (mumol CO2 m-2 s-1)
1484	REAL(r_std) :: Rm !! Day respiration in the mesophyll (umol CO2 m-2 s-1)
1485	REAL(r_std) :: Rdtot !! Total Day respiration (respiratory CO2 release other than
1486	!! by photorespiration) (mumol CO2 m-2 s-1)
1487	REAL(r_std) :: Sco !! Relative CO2 /O2 specificity factor for Rubisco (bar bar-1)
1488	REAL(r_std) :: S_Vcmax_acclim_temp !! Entropy term for Vcmax
1489	REAL(r_std) :: speed !! wind speed @tex ($m s^{-1}$) @endtex
1490	REAL(r_std) :: sum_rveget_rstruct !! recalculated rveget
1491	!! Final unit is @tex ($m s^{-1}$) @endtex
1492	REAL(r_std) :: T_gamma_star !! Temperature dependance of gamma_star (unitless)
1493	REAL(r_std) :: T_gm !! Temperature dependance of gmw
1494	REAL(r_std) :: T_Kmc !! Temperature dependance of KmC (unitless)
1495	REAL(r_std) :: T_KmO !! Temperature dependance of KmO (unitless)
1496	REAL(r_std) :: T_Sco !! Temperature dependance of Sco
1497	REAL(r_std) :: T_Rd !! Temperature dependance of Rd (unitless)
1498	REAL(r_std) :: T_Vcmax !! Temperature dependance of Vcmax (unitless)
1499	REAL(r_std) :: VPD !! Vapor Pressure Deficit (kPa)
1500	!! age (umol CO2 (gN)-1 s-1)
1501	REAL(r_std) :: VpJ2 !! e- transport-limited PEP carboxylation rate (umol CO2 m-2 s-1)
1502	REAL(r_std) :: vcmax !! maximum rate of carboxylation
1503	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
1504	REAL(r_std) :: vc !! Maximum rate of Rubisco activity-limited carboxylation
1505	!! (mumol CO2 mâ2 sâ1)
1506	REAL(r_std) :: vc2 !! rate of carboxylation (at a specific LAI level)
1507	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
1508	REAL(r_std) :: water_lim !! water limitation factor (0-1,unitless)
1509	REAL(r_std) :: zqsvegrap !! relative water quantity in the water
1510	!! interception reservoir (0-1,unitless)
1511	REAL(r_std) :: z !! A lumped parameter (see Yin et al. 2009) ( mol mol-1)
1512	REAL(r_std) :: Cc_loc !! Cellular CO2 concentration
1513	!_ ================================================================================================================================
1514
1515	! 0. Initialization ++ CHECK dimension
1516
1517	lai_total = zero
1518	assimi(:) = zero
1519	gs(:) = zero
1520	gstot = zero
1521	assimtot = zero
1522	Rd(:) = zero
1523	Rdtot=zero
1524	vbeta3pot(ji,jv) = zero
1525	gsmean(ji,jv) = zero
1526	gpp(ji,jv) = zero
1527	light_tran_to_level(:)=un
1528	cimean(ji,jv) = Ca(ji)
1529	vc2 = zero
1530
1531	! 1.1 Photosynthesis parameters
1532	!+++CHECK+++
1533	! Why is total lai calculated as the sum of lai_per_layer rather than the more
1534	! common approach by using ccc_to_lai ? lai(ji,jv) = cc_to_lai( circ_class_biomass(ji,jv,:, ileaf,&
1535	! icarbon),circ_class_n(ji,jv,:),jv)?
1536	! Is there a difference between both approaches? Probably best to test before changing this
1537	! code. I tried to test but the test requires running the multi-layer energy budget which
1538	! crashed in stomate_soilcarbon.f90 see ticket 457.
1539	DO ilev=1,nlevels_tot
1540	!! total LAI for column in question
1541	lai_total = lai_total + lai_per_level(ji,jv,ilev)
1542	ENDDO
1543
1544	! Parametrization of Medlyn et al. (2002) - from Bernacchi et al. (2001)
1545	T_KmC = Arrhenius(temp_air(ji),298.,E_KmC(jv))
1546	T_KmO = Arrhenius(temp_air(ji),298.,E_KmO(jv))
1547	T_Sco = Arrhenius(temp_air(ji),298.,E_Sco(jv))
1548	T_gamma_star = Arrhenius(temp_air(ji),298.,E_gamma_star(jv))
1549	! Parametrization of Yin et al. (2009) - from Bernacchi et al. (2001)
1550	T_Rd = Arrhenius(temp_air(ji),298.,E_Rd(jv))
1551
1552	KmC = KmC25(jv)*T_KmC
1553	KmO = KmO25(jv)*T_KmO
1554	Sco = Sco25(jv)*T_sco
1555	gamma_star = gamma_star25(jv)*T_gamma_star
1556	! low_gamma_star is defined by Yin et al. (2009)
1557	! as the half of the reciprocal of Sco - See Table 2
1558	! We derive its value from Equation 7 of Yin et al. (2009)
1559	! assuming a constant O2 concentration (Oi)
1560	low_gamma_star = gamma_star/Oi
1561
1562	! 1.2 taking N into account for the calculations of Vcmax, as in
1563	! diffuco_trans_co2.
1564	nue = assim_param(ji,jv,inue)
1565	ntot = assim_param(ji,jv,ileafn)
1566	IF (lai_total .GT. min_sechiba) THEN
1567	leafN_top = ntot * ext_coeff_N(jv) / ( 1.-exp(-ext_coeff_N(jv) * lai_total) )
1568	ELSE
1569	leafN_top = zero
1570	ENDIF
1571	vcmax = nue*leafN_top
1572
1573	S_Vcmax_acclim_temp = aSV(jv) + bSV(jv) * MAX(11.,MIN(temp_growth(ji),35.))
1574	T_Vcmax = Arrhenius_modified_nodim(temp_air(ji),298., &
1575	E_Vcmax(jv),D_Vcmax(jv),S_Vcmax_acclim_temp)
1576	vc = vcmax * T_Vcmax
1577
1578	! 1.3 Estimate relative humidity of air
1579	CALL qsatcalc (kjpindex, temp_air, pb, qsatt)
1580
1581	VPD = ( qsatt(ji)pb(ji) / (Tetens_1+qsatt(ji)Tetens_2 ) ) - &
1582	( qair(ji) * pb(ji) / (Tetens_1+qair(ji)* Tetens_2) )
1583	! VPD is needed in kPa
1584	! We limit the impact of VPD in the range of [0:6] kPa
1585	! This seems to be the typical range of values
1586	! that the vapor pressure difference can take, although
1587	! there is a paper by Mott and Parkhurst, Plant, Cell & Environment,
1588	! Vol. 14, pp. 509--519 (1991) which gives a value of 15 kPa.
1589	VPD = MAX(zero,MIN(VPD/10.,6.))
1590	fvpd = 1. / ( 1. / (a1(jv) - b1(jv) * VPD) - 1. )
1591
1592	gb = gb_ref * 44.6 * (tp_00/temp_air(ji)) * (pb(ji)/pb_std)
1593	! conversion from (mol H2O m-2 s-1) to (mol CO2 m-2 s-1)
1594	gb = gb / ratio_H2O_to_CO2
1595
1596	! 1.4 Intermediate variables needed for recalculation
1597	jwind = SQRT (u(ji)u(ji) + v(ji)v(ji))
1598	speed = MAX(min_wind, jwind)
1599
1600	zqsvegrap = zero
1601	IF (qsintmax(ji,jv) .GT. min_sechiba) THEN
1602	!! relative water quantity in the water interception reservoir
1603	zqsvegrap = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv))
1604	ENDIF
1605
1606	T_gm = Arrhenius_modified_nodim(temp_air(ji),298.,E_gm(jv),D_gm(jv),S_gm(jv))
1607	! water_lim is always 1 for this module, because it in ok_hydrol_arch already
1608	water_lim = un
1609	gm = gm25(jv) * T_gm * MAX(1-stress_gm(jv), water_lim)
1610	g0var = g0(jv)* MAX(1-stress_gs(jv), water_lim)
1611
1612	! give a default value of ci for all pixel that do not assimilate
1613	DO ilev=1,nlevels_tot
1614	new_leaf_ci(ji,jv,ilev) = Ca(ji)
1615	ENDDO
1616
1617	! 2. Recalculating assimilation
1618
1619	! 2.1 Recalculating resistance
1620	! recalculation_1) creist
1621	! If the leaf has a full film of water then the original equation for
1622	! vbeta3 doesn't include cresist in it so the connection is lost.
1623	! In that case, cresist is recalculated using a simple ratio of a reduction of
1624	! vbeta3. Note that this approach follows the way of calculation in diffuco but
1625	! overlooks stomatal closure through water film.
1626	IF (zqsvegrap .EQ. un) THEN
1627	cresist(ji,jv) = cresist(ji,jv) *ratio_vbeta3(ji,jv)
1628	ELSE
1629	cresist(ji,jv) = MAX(vbeta3(ji,jv) /veget(ji,jv) , &
1630	(vbeta3(ji,jv) - vbeta23(ji,jv)) /(un - zqsvegrap)/veget(ji,jv))
1631	ENDIF
1632
1633	! recalculation_2) sum_rveget_rstruct
1634	IF (cresist(ji,jv) .GT. zero .AND. speed .GT. min_sechiba .AND. &
1635	q_cdrag(ji) .GT. min_sechiba .AND. rveg_pft(jv) .GT. min_sechiba) THEN
1636
1637	sum_rveget_rstruct = (un/cresist(ji,jv) - un) / (speed * q_cdrag(ji) * &
1638	(veget(ji,jv)/veget_max(ji,jv)) * rveg_pft(jv))
1639	ELSE
1640	sum_rveget_rstruct = 9999.
1641	END IF
1642
1643	vbeta3pot(ji,jv) = MAX(zero, veget_max(ji,jv) *cresist(ji,jv))
1644
1645	! recalculation_3) gstot
1646	IF (sum_rveget_rstruct .ne. 9999.) THEN
1647	gstot = un / sum_rveget_rstruct
1648	ELSE
1649	! Need to convert the units of g0
1650	gstot = g0(jv) * mol_to_m_1 * (temp_air(ji)/tp_00) * &
1651	(pb_std/pb(ji))*ratio_H2O_to_CO2
1652	END IF
1653
1654	! One equation normliznig gstot using
1655	! lai_coupled is omitted because lai_coupled is the same with lai at this
1656	! moment. It should be applied if the change in lai_coupled is made.
1657
1658	! we now do a unit conversion, from m/s to mol/m^2/s, the reverse
1659	! of the calculation of section 2.5,
1660
1661	gstot = gstot / (mol_to_m_1 (temp_air(ji)/tp_00) &
1662	(pb_std/pb(ji))*ratio_H2O_to_CO2)
1663
1664	DO ilev = nlevels_tot-1,1,-1
1665	light_tran_to_level(ilev)=light_tran_to_level(ilev+1)*Light_Tran_Tot(ji,jv,ilev+1)
1666	ENDDO
1667
1668	DO ilev = 1, nlevels_tot
1669	IF (Light_Abs_Tot(ji,jv,ilev) .GT. min_sechiba) THEN
1670	IF(lai_per_level(ji,jv,ilev) .LE. min_sechiba)THEN
1671	! This is not necessarily a problem. LAI is updated in slowproc, whil
1672	! absorbed light is updated in condveg, which is called after diffuco.
1673	! So we might be one timestep off. It generally should be at night,
1674	! though, so there is no absorbed light. If this persists, that
1675	! is a real problem.
1676	IF(printlev_loc >= 3)THEN
1677	WRITE(numout,*) 'WARNING: We have absorbed light but no LAI ' // &
1678	'in this level. Skipping.'
1679	WRITE(numout,'(A,3I8)') 'WARNING: ji,jv,ilev: ',ji,jv,ilev
1680	WRITE(numout,'(A,2F20.10)') 'WARNING: Light_Abs_Tot(ji,jv,ilev), ' // &
1681	'lai_per_level(ji,jv,ilev): ',&
1682	Light_Abs_Tot(ji,jv,ilev), lai_per_level(ji,jv,ilev)
1683	ENDIF
1684	CYCLE
1685	ENDIF ! lai_per_level
1686
1687	!
1688	! nitrogen is scaled according to the transmitted light
1689	! before the caclulation was (1 - exp(lai)) we replaced
1690	! it by the transmitted light calculated by the 2stream
1691	! model, maybe ( un - Light_Tran_Tot(iainia,jv,ilev) )
1692	! could be replaced by the directly calculated absorbed
1693	! light Light_Abs_Tot(iainia,jv,ilev).
1694	! Notice that light in the old code was the cumulative light transmitted to
1695	! the layer, while Light_Tran_Tot is the relative light transmitted
1696	! through the layer. Therefore, we use a different variable.
1697	!
1698	N_profile = ( un - .7_r_std * ( un - light_tran_to_level(ilev) ) )
1699	! see Comment in legend of Fig. 6 of Yin et al. (2009)
1700	! Rd25 is assumed to equal 0.01 Vcmax25
1701	! Effects of drought stress on respiration are less
1702	! well known: increases, decreases and no change have been found.
1703	! For now there is no effect of drought stress on dark respiration
1704	! when hydrol_arch is used.
1705	vc2 = vc * N_profile * MAX(1-stress_vcmax(jv), water_lim)
1706	Rd(ilev) = vcmax * N_profile * 0.01 * T_Rd
1707
1708	IF (printlev_loc >= 4) THEN
1709	WRITE(numout,*) 'mleb_diffuco_trans_co2: recalculation starts'
1710	ENDIF
1711
1712	! 2.2 Assimilation for C4 plants
1713	IF ( is_c4(jv) ) THEN
1714
1715	! Analytical resolution of the Assimilation based
1716	! Yin et al. (2009)
1717	! Eq. 28 of Yin et al. (2009)
1718	fcyc= 1. - ( 4.(1.-fpsir(jv))(1.+fQ(jv)) + &
1719	3.h_protons(jv)fpseudo(jv) ) / &
1720	( 3.h_protons(jv) - 4.(1.-fpsir(jv)))
1721
1722	! See paragraph after eq. (20b) of Yin et al.
1723	Rm=Rd(ilev)/2.
1724
1725	! We assume that cs_star equals ci_star
1726	! (see Comment in legend of Fig. 6 of Yin et al. (2009)
1727	! Equation 26 of Yin et al. (2009)
1728	Cs_star = (gbs(jv) * low_gamma_star * Oi - ( 1. + &
1729	low_gamma_star * alpha(jv) / 0.047) * Rd(ilev) + Rm )/ &
1730	( gbs(jv) + kp(jv) )
1731
1732	! eq. 11 of Yin et al (2009)
1733	J2 = JJ_out(ji,jv,ilev) / ( 1. - fpseudo(jv) / ( 1. - fcyc ) )
1734
1735	! Equation right after eq. (20d) of Yin et al. (2009)
1736	z = ( 2. + fQ(jv) - fcyc ) / ( h_protons(jv) * (1. - fcyc ))
1737
1738	VpJ2 = fpsir(jv) * J2 * z / 2.
1739	! See paragraph after eq. (20b) of Yin et al.
1740	Rm=Rd(ilev)/2.
1741
1742	! We assume that cs_star equals ci_star (see Comment in legend of Fig. 6 of Yin et al. (2009)
1743	! Equation 26 of Yin et al. (2009)
1744	IF ( info_limitphoto(ji,jv,ilev) .EQ. 1 ) THEN
1745	assim_coeff_a = 1. + kp(jv) / gbs(jv)
1746	assim_coeff_b = 0.
1747	assim_coeff_x1 = vc2
1748	assim_coeff_x2 = KmC/KmO
1749	assim_coeff_x3 = KmC
1750	! Is J limiting the Assimilation
1751	ELSE IF (info_limitphoto(ji,jv,ilev) .EQ. 2) THEN
1752	assim_coeff_a = 1.
1753	assim_coeff_b = VpJ2
1754	assim_coeff_x1 = (1.- fpsir(jv)) * J2 * z / 3.
1755	assim_coeff_x2 = 7. * low_gamma_star / 3.
1756	assim_coeff_x3 = 0.
1757	ELSE
1758	IF ( printlev_loc >= 4) THEN
1759	WRITE(numout,*) 'limit_photo was not defined: no assimilation happened'
1760	WRITE(numout,*) 'ji,jv,ilev',ji,jv,ilev
1761	ENDIF
1762	assimi(ilev) = -Rd(ilev)
1763	CYCLE
1764	ENDIF
1765
1766	! recalculation_4) gs
1767	gs(ilev) = gstot*gstot_frac(ji,jv,ilev) / lai_per_level(ji,jv,ilev)
1768
1769	IF ( gs(ilev) .LT. min_sechiba ) THEN
1770	gs(ilev) = g0(jv)
1771	ENDIF
1772	! IF (gs(iainia) .gt. gstot(iainia)) THEN
1773	! gs(iainia) = gstot(iainia)
1774	! END IF
1775
1776	IF (gs(ilev) .GE. g0(jv)) THEN
1777	! recalculation_5) assimi (C4) Equations copied from
1778	! diffuco_trans_co2 in diffuco.f90 and resuffled to isolate
1779	! assimi
1780	assimi(ilev) = ( (gs(ilev)-g0(jv))(Ca(ji)-Cs_star) - Rd(ilev)fvpd ) / &
1781	( fvpd +(gs(ilev)-g0(jv))/gb )
1782	ELSE
1783
1784	assimi(ilev) = -Rd(ilev)
1785	ENDIF
1786
1787	! recalculation_6) Ci. Equations copied from diffuco_trans_co2 in diffuco.f90
1788	Obs = ( alpha(jv) * assimi_lev(ji,jv,ilev) ) / ( 0.047 * gbs(jv) ) + Oi
1789	! Eq. 23 of Yin et al. (2009)
1790	Cc_loc = (( assimi_lev(ji,jv,ilev) + Rd(ilev)) * &
1791	( assim_coeff_x2 * Obs + assim_coeff_x3 ) +&
1792	low_gamma_star * Obs * assim_coeff_x1 )/ &
1793	MAX(min_sechiba, assim_coeff_x1 - ( assimi_lev(ji,jv,ilev) + Rd(ilev) ))
1794	! Eq. 22 of Yin et al. (2009)
1795	new_leaf_ci(ji,jv,ilev) = ( Cc_loc - ( assim_coeff_b - assimi(ilev) - Rm ) / &
1796	gbs(jv) ) / assim_coeff_a
1797
1798	IF(printlev_loc >=4) THEN
1799	WRITE(numout,*) 'after recalculation'
1800	WRITE(numout,*) 'ji,jv,ilev',ji,jv,ilev
1801	WRITE(numout,*) 'cresist',cresist(ji,jv)
1802	WRITE(numout,*) 'sum_rveget_rstruct',sum_rveget_rstruct
1803	WRITE(numout,*) 'gs',gs(ilev)
1804	WRITE(numout,*) 'assimi',assimi(ilev)
1805	WRITE(numout,*) 'Rd', Rd(ilev)
1806	ENDIF
1807
1808	IF (assimi(ilev) .lt. zero) THEN
1809	assimi(ilev) = zero
1810	! WRITE(numout,) 'assimi(iainia,jv) is negative **************************************'
1811	! WRITE(numout,*) 'assimi(iainia,jv) is: ', assimi(iainia,jv)
1812	! WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
1813	! WRITE(numout,*) 'g0(jv) is: ', g0(jv)
1814	! WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
1815	! WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
1816	! WRITE(numout,*) 'Cs_star is: ', Cs_star
1817	! WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
1818	! WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
1819	! WRITE(numout,) '********************************************************'
1820	ELSE IF (assimi(ilev) .gt. 500.) THEN
1821	! WRITE(numout,*) 'assimi(iainia,jv) is very large'
1822	! WRITE(numout,*) 'assimi(iainia,jv) is: ', assimi(iainia,jv)
1823	! WRITE(numout,*) 'temp_air(iainia) is: ', temp_air(iainia)
1824	! WRITE(numout,*) 'vbeta3(iainia, jv) is: ', vbeta3(iainia, jv)
1825	! WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
1826	! WRITE(numout,*) 'g0(jv) is: ', g0(jv)
1827	! WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
1828	! WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
1829	! WRITE(numout,*) 'Cs_star is: ', Cs_star
1830	! WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
1831	! WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
1832	! WRITE(numout,*) 'gstot(iainia) is: ', gstot(iainia)
1833	! WRITE(numout,*) 'sum_rveget_rstruct(iainia,jv) is: ', sum_rveget_rstruct(iainia,jv)
1834	! WRITE(numout,*) 'speed is: ', speed
1835	! WRITE(numout,*) 'q_cdrag(iainia) is: ',q_cdrag(iainia)
1836	! WRITE(numout,*) 'rveg_pft(jv) is: ', rveg_pft(jv)
1837	! WRITE(numout,*) 'gstot_frac(iainia,jv,ilev) is: ', gstot_frac(iainia,jv,ilev)
1838	! WRITE(numout,*) 'lai_per_level(iainia,jv,ilev) is: ', lai_per_level(iainia,jv,ilev)
1839	! WRITE(numout,) '********************************************************'
1840	END IF
1841
1842
1843	ELSE !is_c4
1844
1845	! 2.2 Assimilation for C3 plants
1846	! Eq. 17 of Yin et al. (2009)
1847	! See eq. right after eq. 15 of Yin et al. (2009)
1848	ci_star = gamma_star
1849
1850	IF ( info_limitphoto(ji,jv,ilev) .EQ. 1 ) THEN
1851	assim_coeff_x1 = vc2
1852	! It should be O not Oi (comment from Vuichard)
1853	assim_coeff_x2 = KmC * ( 1. + 2gamma_starSco/ KmO )
1854
1855	! Is J limiting the Assimilation
1856	ELSE IF ( info_limitphoto(ji,jv,ilev) .EQ. 2) THEN
1857	assim_coeff_x1 = JJ_out(ji,jv,ilev)/4.
1858	assim_coeff_x2 = 2. * gamma_star
1859
1860	ELSE
1861	IF ( printlev_loc >= 4) THEN
1862	WRITE(numout,*) 'limit_photo was not defined: no assimilation happened'
1863	WRITE(numout,*) 'ji,jv,ilev',ji,jv,ilev
1864	ENDIF
1865	assimi(ilev) = -Rd(ilev)
1866	CYCLE
1867	ENDIF
1868
1869	! recalculation_4) gs
1870	gs(ilev) = gstot * gstot_frac(ji,jv,ilev) / lai_per_level(ji,jv,ilev)
1871	! IF (gs(iainia) .gt. gstot(iainia)) THEN
1872	! gs(iainia) = gstot(iainia)
1873	! END IF
1874	IF ( gs(ilev) .LT. min_sechiba ) THEN
1875	gs(ilev) = g0(jv)
1876	ENDIF
1877
1878
1879	IF(gs(ilev) .GE. g0(jv)) THEN
1880
1881	! recalcaultion_5) assimi (C3)
1882	! Unlike assimi for c4 plants, assimi calculation for c3 plants contains Cc_loc
1883	! inside which needs 5 equations to combine including the cubic equation
1884	! (Yin and Struik (2009), Eq 15 to 19), which can, in turn, make the
1885	! recalculation quite complex. To avoid that, only Eq 15 and 16 were used
1886	! assuming that gm and Cc_loc will be constant. The result from the simplified
1887	! (below) equation is not 100% precise but precise enough to accept
1888	! (~0.3% lower) when it is checked in diffuco.
1889	!
1890	assimi(ilev) =(gb * gs(ilev) * (Ca(ji) * (gs(ilev) - g0(jv)) + &
1891	gamma_star * (g0(jv) - gs(ilev)) -Rd(ilev)*fvpd))/ &
1892	(gb * ( -g0(jv) +gs(ilev)*fvpd + gs(ilev)) + &
1893	gs(ilev)*(gs(ilev) - g0(jv)))
1894	ELSE
1895	assimi(ilev) = -Rd(ilev)
1896	IF(printlev_loc >=4) THEN
1897	WRITE(numout,*) 'gs < g0. assimi will be set to -Rd'
1898	WRITE(numout,*) 'ji,jv,ilev',ji,jv,ilev
1899	ENDIF
1900	ENDIF
1901
1902
1903	! recalculation_6) Ci
1904	! Recalculated Ci is too low!! Decrease in assimi decreases Cc_loc way too much, which
1905	! results in too low leaf_ci. It can give a huge impact on C13, so
1906	! C13=y is prohibited for now. If assimi_lev is used instead of
1907	! assimi for Cc_loc calculation (keeping the Cc_loc without stress) the change in
1908	! leaf_ci is reasonable, but it was not applied because ci is not
1909	! written in restart file and it doesn't seem right.
1910
1911	Cc_loc = (gamma_star * assim_coeff_x1 + (assimi(ilev)+Rd(ilev)) * assim_coeff_x2) / &
1912	max(min_sechiba, assim_coeff_x1 - (assimi(ilev)+Rd(ilev)))
1913	new_leaf_ci(ji,jv,ilev) = Cc_loc + assimi(ilev)/gm
1914
1915
1916	IF(test_pft == jv .AND. printlev_loc >=4) THEN
1917	WRITE(numout,*) 'after recalculation, C3'
1918	WRITE(numout,*) 'ji,jv,ilev',ji,jv,ilev
1919	WRITE(numout,*) 'cresist',cresist(ji,jv)
1920	WRITE(numout,*) 'sum_rveget_rstruct',sum_rveget_rstruct
1921	WRITE(numout,*) 'gs',gs(ilev)
1922	WRITE(numout,*) 'assimi_recal',assimi(ilev)
1923	WRITE(numout,*) 'Rd', Rd(ilev)
1924	WRITE(numout,*) 'Cc_recal',Cc_loc
1925	WRITE(numout,*) 'Ci_recal',new_leaf_ci(ji,jv,ilev)
1926	ENDIF
1927
1928	IF (assimi(ilev) .lt. min_sechiba) THEN
1929	assimi(ilev) = zero
1930	! WRITE(numout,) 'assimi(iainia,jv) is negative **************************************'
1931	! WRITE(numout,*) 'assimi(iainia,jv) is: ', assimi(iainia,jv)
1932	! WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
1933	! WRITE(numout,*) 'g0(jv) is: ', g0(jv)
1934	! WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
1935	! WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
1936	! WRITE(numout,*) 'Cs_star is: ', Cs_star
1937	! WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
1938	! WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
1939	! WRITE(numout,) '********************************************************'
1940	ELSE IF (assimi(ilev) .gt. 500.) THEN
1941	! WRITE(numout,*) 'assimi(iainia,jv) is very large'
1942	! WRITE(numout,*) 'assimi(iainia,jv) is: ', assimi(iainia,jv)
1943	! WRITE(numout,*) 'temp_air(iainia) is: ', temp_air(iainia)
1944	! WRITE(numout,*) 'vbeta3(iainia, jv) is: ', vbeta3(iainia, jv)
1945	! WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
1946	! WRITE(numout,*) 'g0(jv) is: ', g0(jv)
1947	! WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
1948	! WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
1949	! WRITE(numout,*) 'Cs_star is: ', Cs_star
1950	! WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
1951	! WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
1952	! WRITE(numout,*) 'gstot(iainia) is: ', gstot(iainia)
1953	! WRITE(numout,*) 'sum_rveget_rstruct(iainia,jv) is: ', sum_rveget_rstruct(iainia,jv)
1954	! WRITE(numout,*) 'speed is: ', speed
1955	! WRITE(numout,*) 'q_cdrag(iainia) is: ',q_cdrag(iainia)
1956	! WRITE(numout,*) 'rveg_pft(jv) is: ', rveg_pft(jv)
1957	! WRITE(numout,*) 'gstot_frac(iainia,jv,ilev) is: ', gstot_frac(iainia,jv,ilev)
1958	! WRITE(numout,*) 'lai_per_level(iainia,jv,ilev) is: ', lai_per_level(iainia,jv,ilev)
1959	! WRITE(numout,) '********************************************************'
1960	END IF
1961	ENDIF !is_c4
1962
1963	IF(test_pft == jv .AND. printlev_loc>=4)THEN
1964	WRITE(numout,*) 'assimi ',ilev,ji,assimi(ilev)
1965	WRITE(numout,*) 'lai_per_level',lai_per_level(ji,jv,ilev)
1966	ENDIF
1967
1968	ENDIF !IF Light
1969	ENDDO !DO ilev
1970
1971	!Initialize gstot to calculate again with recalculated gs
1972	gstot = zero
1973	!! 2.3 Integration at the canopy level
1974	DO ilev = 1,nlevels_tot
1975	assimtot = assimtot + assimi(ilev) * lai_per_level(ji,jv,ilev)
1976	Rdtot = Rdtot + &
1977	Rd(ilev) * lai_per_level(ji,jv,ilev)
1978	gstot = gstot + &
1979	gs(ilev) * lai_per_level(ji,jv,ilev)
1980	ENDDO
1981
1982	IF (err_act.GT.1 .AND. assimtot .gt. 20) THEN
1983	WRITE(numout,*) 'WARNING: assimtot is very high'
1984	WRITE(numout,*) 'assimtot is: ', assimtot
1985	DO ilev=1,nlevels_tot
1986	WRITE (numout,*) 'level: ', ilev, ' pft: ', jv , ' assimi_store: ', assimi(ilev), &
1987	&' lai_per_level: ', lai_per_level(ji, jv, ilev)
1988	END DO
1989	END IF
1990
1991	!! 2.4 Calculate resistances
1992	! Copy from diffuco.f90
1993
1994	!! Mean stomatal conductance for CO2 (mol m-2 s-1)
1995	gsmean(ji,jv) = gstot
1996
1997	! cimean is completely diagnostic and not
1998	! used for anything. If the analytical photosynthesis is unable to find
1999	! a solution for the A_1 coefficients in any of the levels, assimtot will
2000	! be equal to -Rd, which means we will be dividing by zero. This loop
2001	! is to prevent that, and since cimean is diagnostic is doesn't matter
2002	! what the value is.
2003	IF ( ABS(gsmean(ji,jv)-g0var*lai_total) .GT. min_sechiba) THEN
2004	cimean(ji,jv) = (fvpd*(assimtot+Rdtot)) /&
2005	(gsmean(ji,jv)-g0var*lai_total) + gamma_star
2006	ELSE
2007	cimean(ji,jv) = gamma_star
2008	ENDIF
2009
2010	gpp(ji,jv) = assimtot12e-6veget_max(ji,jv)*dt_sechiba
2011
2012	END SUBROUTINE mleb_diffuco_trans_co2_short
2013
2014
2015	!! ADD HEADER
2016
2017	!_ =================================================================================================================================
2018
2019	SUBROUTINE jnetcdf_trans_co2(kjpindex, james_netcdf_flag, gs_distribution_in, gs_output_in, &
2020	lai_per_level_in, light_tran_to_level_in, gs_diffuco_output_in, &
2021	gstot_output_in)
2022
2023	USE netcdf
2024
2025	IMPLICIT NONE
2026
2027	!! 0.1 Input variables
2028
2029	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
2030
2031	CHARACTER(LEN=80) :: jamescdfname
2032
2033	REAL(r_std), DIMENSION (nvm, nlevels_tot), INTENT(in) :: gs_distribution_in
2034	REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in) :: gs_output_in
2035	REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in) :: gstot_output_in
2036	REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in) :: lai_per_level_in
2037	REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in) :: light_tran_to_level_in
2038	REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in) :: gs_diffuco_output_in
2039
2040
2041	LOGICAL :: james_netcdf_flag
2042
2043
2044	!! 0.2 Output variables
2045
2046	!! 0.3 Modified variables
2047
2048	!! 0.4 Local variables
2049
2050	INTEGER(i_std), SAVE :: gsd_id
2051	!$OMP THREADPRIVATE(gsd_id)
2052	INTEGER(i_std), SAVE :: gso_id
2053	!$OMP THREADPRIVATE(gso_id)
2054	INTEGER(i_std), SAVE :: gsto_id
2055	!$OMP THREADPRIVATE(gsto_id)
2056	INTEGER(i_std), SAVE :: lpl_id
2057	!$OMP THREADPRIVATE(lpl_id)
2058	INTEGER(i_std), SAVE :: lttl_id
2059	!$OMP THREADPRIVATE(lttl_id)
2060	INTEGER(i_std), SAVE :: gsdo_id
2061	!$OMP THREADPRIVATE(gsdo_id)
2062
2063	INTEGER(i_std) :: ii, jj
2064	INTEGER(i_std), SAVE :: mleb_step_count
2065	!$OMP THREADPRIVATE(mleb_step_count)
2066
2067	INTEGER(i_std), SAVE :: iret, iret2, fid, fid2
2068	!$OMP THREADPRIVATE(iret,iret2,fid,fid2)
2069	INTEGER(i_std) :: TimeDimID, jtid, ts_id, HeightDimID, i
2070	INTEGER(i_std) :: height_dim_id, height_dim_id2, time_dim_id
2071	INTEGER(i_std) :: pftno_dim_id
2072	INTEGER(i_std) :: dimids(3)
2073	INTEGER(i_std), SAVE :: tap_id
2074	!$OMP THREADPRIVATE(tap_id)
2075	INTEGER, DIMENSION(3) :: jstartnc, jcountnc
2076
2077
2078	jamescdfname = '/home/users/jryder/james.ryder/outputcdf.nc'
2079
2080	IF (james_netcdf_flag) THEN
2081
2082	james_netcdf_flag = .FALSE.
2083
2084	! mleb_step_count = 0
2085
2086	! WRITE(numout,*) 'james_netcdf_flag call here'
2087
2088	! âŠ define the parameters here
2089
2090	iret = nf90_create(jamescdfname, nf90_write, fid)
2091
2092
2093	! define the dimensions
2094
2095	iret = nf90_def_dim(fid, "pftno", nvm, pftno_dim_id)
2096	iret = nf90_def_dim(fid, "height", nlevels_tot, height_dim_id)
2097	iret = nf90_def_dim(fid, "time", nf90_unlimited, time_dim_id)
2098
2099	dimids = (/ pftno_dim_id, height_dim_id, time_dim_id /)
2100
2101
2102	! define the variables
2103	iret = nf90_def_var (fid, "gs_distribution_in", nf90_real4, dimids, gsd_id)
2104	iret = nf90_def_var (fid, "gs_output_in", nf90_real4, dimids, gso_id)
2105	iret = nf90_def_var (fid, "gstot_output_in", nf90_real4, dimids, gsto_id)
2106	iret = nf90_def_var (fid, "lai_per_level_in", nf90_real4, dimids, lpl_id)
2107	iret = nf90_def_var (fid, "light_tran_to_level_in", nf90_real4, dimids, lttl_id)
2108	iret = nf90_def_var (fid, "gs_diffuco_output_in", nf90_real4, dimids, gsdo_id)
2109
2110
2111	iret = nf90_enddef (fid)
2112
2113	END IF !(james_netcdf_flag)
2114
2115
2116	! open the netCDF file at each step
2117	iret = nf90_open(jamescdfname, nf90_Write, fid)
2118
2119	jstartnc = (/ 1, 1, 1 /) ! vector of integers specifying the index in the variable from
2120	! which the first of the data values will be read
2121	jcountnc = (/ nvm, nlevels_tot, 1 /) ! vector of integers specifying the number of indices
2122	! selected along each dimension
2123
2124	mleb_step_count = mleb_step_count + 1 ! this is the time step counter
2125
2126
2127	jstartnc(3) = mleb_step_count ! (2) refers to the time dimension at present (it
2128	! must always come last)
2129
2130
2131	iret = nf90_put_var(fid, gsd_id, gs_distribution_in, start = jstartnc, count = jcountnc)
2132	iret = nf90_put_var(fid, gso_id, gs_output_in, start = jstartnc, count = jcountnc)
2133	iret = nf90_put_var(fid, gsto_id, gstot_output_in, start = jstartnc, count = jcountnc)
2134	iret = nf90_put_var(fid, lpl_id, lai_per_level_in, start = jstartnc, count = jcountnc)
2135	iret = nf90_put_var(fid, lttl_id, light_tran_to_level_in, start = jstartnc, count = jcountnc)
2136	iret = nf90_put_var(fid, gsdo_id, gs_diffuco_output_in, start = jstartnc, count = jcountnc)
2137
2138
2139
2140	! close the netCDF file and write to disk
2141	iret = nf90_close(fid)
2142
2143
2144	END SUBROUTINE jnetcdf_trans_co2
2145
2146
2147	! -------------------------------------------------------------------------------------------------------------------------------------
2148
2149
2150
2151	END MODULE hydraulic_arch

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_sechiba/hydraulic_arch.f90 @ 7852

Download in other formats: