Context Navigation

source: branches/publications/ORCHIDEE_CN_CAN_r5698/src_sechiba/diffuco.f90 @ 7346

Last change on this file since 7346 was 5656, checked in by sebastiaan.luyssaert, 6 years ago
DEV: tested for a single pixel with 5 different configurations to start and 5 different configurations to restart the model. Passed all basic tests for impose_cn = y. These changes contribute to ticket #286 and contain major cleaning and restructuring of slowproc.f90
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 185.9 KB

Line
1	! ================================================================================================================================
2	! MODULE : diffuco
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 limiting coefficients, both aerodynamic
10	!! and hydrological, for the turbulent heat fluxes.
11	!!
12	!!\n DESCRIPTION: The aerodynamic resistance R_a is used to limit
13	!! the transport of fluxes from the surface layer of vegetation to the point in the atmosphere at which
14	!! interaction with the LMDZ atmospheric circulation model takes place. The aerodynamic resistance is
15	!! calculated either within the module r_aerod (if the surface drag coefficient is provided by the LMDZ, and
16	!! if the flag 'ldq_cdrag_from_gcm' is set to TRUE) or r_aero (if the surface drag coefficient must be calculated).\n
17	!!
18	!! Within ORCHIDEE, evapotranspiration is a function of the Evaporation Potential, but is modulated by a
19	!! series of resistances (canopy and aerodynamic) of the surface layer, here represented by beta.\n
20	!!
21	!! DESCRIPTION :
22	!! \latexonly
23	!! \input{diffuco_intro.tex}
24	!! \endlatexonly
25	!! \n
26	!!
27	!! This module calculates the beta for several different scenarios:
28	!! - diffuco_snow calculates the beta coefficient for sublimation by snow,
29	!! - diffuco_inter calculates the beta coefficient for interception loss by each type of vegetation,
30	!! - diffuco_bare calculates the beta coefficient for bare soil,
31	!! - diffuco_trans_co2 calculates the beta coefficient for transpiration for each type of vegetation, using Farqhar's formula
32	!! - chemistry_bvoc calculates the beta coefficient for emissions of biogenic compounds \n
33	!!
34	!! Finally, the module diffuco_comb computes the combined $\alpha$ and $\beta$ coefficients for use
35	!! elsewhere in the module. \n
36
37	!! RECENT CHANGE(S): Nathalie le 28 mars 2006 - sur proposition de Fred Hourdin, ajout
38	!! d'un potentiometre pour regler la resistance de la vegetation (rveg is now in pft_parameters)
39	!! October 2018: Removed diffuco_trans using Jarvis formula for calculation of beta coefficient
40	!!
41	!! REFERENCE(S) : None
42	!!
43	!! SVN :
44	!! $HeadURL$
45	!! $Date$
46	!! $Revision$
47	!! \n
48	!_ ================================================================================================================================
49
50	MODULE diffuco
51
52	! modules used :
53	USE function_library, ONLY : cc_to_lai, biomass_to_lai, biomass_to_coupled_lai
54	USE constantes
55	USE constantes_soil
56	USE qsat_moisture
57	USE sechiba_io_p
58	USE ioipsl
59	USE pft_parameters
60	USE grid
61	USE ioipsl_para
62	USE xios_orchidee
63	USE chemistry, ONLY : chemistry_initialize, chemistry_bvoc, chemistry_clear
64	IMPLICIT NONE
65
66	! public routines :
67	PRIVATE
68	PUBLIC :: diffuco_main, diffuco_initialize, diffuco_finalize, diffuco_clear
69
70	INTERFACE Arrhenius_modified
71	MODULE PROCEDURE Arrhenius_modified_0d, Arrhenius_modified_1d
72	END INTERFACE
73
74	!
75	! variables used inside diffuco module : declaration and initialisation
76	!
77	!INTEGER(i_std), SAVE :: printlev_loc !! local printlev for this module
78	!!$OMP THREADPRIVATE(printlev_loc)
79	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: leaf_ci !! intercellular CO2 concentration (ppm)
80	!$OMP THREADPRIVATE(leaf_ci)
81	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: wind !! Wind module (m s^{-1})
82	!$OMP THREADPRIVATE(wind)
83
84	CONTAINS
85
86
87	!! =============================================================================================================================
88	!! SUBROUTINE: diffuco_initialize
89	!!
90	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
91	!!
92	!! DESCRIPTION: Allocate module variables, read from restart file or initialize with default values.
93	!! Call chemistry_initialize for initialization of variables needed for the calculations of BVOCs.
94	!!
95	!! RECENT CHANGE(S): None
96	!!
97	!! REFERENCE(S): None
98	!!
99	!! FLOWCHART: None
100	!! \n
101	!_ ==============================================================================================================================
102	SUBROUTINE diffuco_initialize (kjit, kjpindex, index, &
103	rest_id, lalo, neighbours, resolution, &
104	rstruct, q_cdrag)
105
106	!! 0. Variable and parameter declaration
107	!! 0.1 Input variables
108	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
109	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
110	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-)
111	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (-)
112	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
113	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb),INTENT (in):: neighbours !! Vector of neighbours for each
114	REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
115
116	!! 0.2 Output variables
117	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! Structural resistance for the vegetation
118
119	!! 0.3 Modified variables
120	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
121
122	!! 0.4 Local variables
123	INTEGER :: ilai
124	INTEGER :: jv
125	INTEGER :: ier
126	CHARACTER(LEN=4) :: laistring
127	CHARACTER(LEN=80) :: var_name
128	!_ ================================================================================================================================
129
130	!! Initialize local printlev
131	!printlev_loc=get_printlev('diffuco')
132
133	!! 1. Define flag ldq_cdrag_from_gcm. This flag determines if the cdrag should be taken from the GCM or be calculated.
134	!! The default value is true if the q_cdrag variables was already initialized. This is the case when coupled to the LMDZ.
135
136	!Config Key = CDRAG_FROM_GCM
137	!Config Desc = Keep cdrag coefficient from gcm.
138	!Config If = OK_SECHIBA
139	!Config Def = y
140	!Config Help = Set to .TRUE. if you want q_cdrag coming from GCM (if q_cdrag on initialization is non zero).
141	!Config Keep cdrag coefficient from gcm for latent and sensible heat fluxes.
142	!Config Units = [FLAG]
143	IF ( ABS(MAXVAL(q_cdrag)) .LE. EPSILON(q_cdrag)) THEN
144	ldq_cdrag_from_gcm = .FALSE.
145	ELSE
146	ldq_cdrag_from_gcm = .TRUE.
147	ENDIF
148	CALL getin_p('CDRAG_from_GCM', ldq_cdrag_from_gcm)
149	IF (printlev_loc>=2) WRITE(numout,*) "ldq_cdrag_from_gcm = ",ldq_cdrag_from_gcm
150
151	!! 2. Allocate module variables
152	ALLOCATE (leaf_ci(kjpindex,nvm,nlai),stat=ier)
153	IF (ier /= 0) CALL ipslerr_p(3,'diffuco_initialize','Problem in allocate of variable leaf_ci','','')
154
155	ALLOCATE (wind(kjpindex),stat=ier)
156	IF (ier /= 0) CALL ipslerr_p(3,'diffuco_initialize','Problem in allocate of variable wind','','')
157
158	!! 3. Read variables from restart file
159	IF (printlev>=3) WRITE (numout,*) 'Read DIFFUCO variables from restart file'
160
161	CALL ioconf_setatt_p('UNITS', 's/m')
162	CALL ioconf_setatt_p('LONG_NAME','Structural resistance')
163	CALL restget_p (rest_id, 'rstruct', nbp_glo, nvm, 1, kjit, .TRUE., rstruct, "gather", nbp_glo, index_g)
164	IF ( ALL(rstruct(:,:) == val_exp) ) THEN
165	DO jv = 1, nvm
166	rstruct(:,jv) = rstruct_const(jv)
167	ENDDO
168	ENDIF
169
170	CALL ioconf_setatt_p('UNITS', 'ppm')
171	CALL ioconf_setatt_p('LONG_NAME','Leaf CO2')
172	CALL restget_p (rest_id, 'leaf_ci', nbp_glo, nvm, nlai, kjit, .TRUE.,leaf_ci, "gather", nbp_glo, index_g)
173
174	!Config Key = DIFFUCO_LEAFCI
175	!Config Desc = Initial leaf CO2 level if not found in restart
176	!Config If = OK_SECHIBA
177	!Config Def = 233.
178	!Config Help = The initial value of leaf_ci if its value is not found
179	!Config in the restart file. This should only be used if the model is
180	!Config started without a restart file.
181	!Config Units = [ppm]
182	CALL setvar_p (leaf_ci, val_exp,'DIFFUCO_LEAFCI', 233._r_std)
183
184
185	!! 4. Initialize chemistry module
186	IF (printlev>=3) WRITE(numout,*) "ok_bvoc:",ok_bvoc
187	IF ( ok_bvoc ) CALL chemistry_initialize(kjpindex, lalo, neighbours, resolution)
188
189	END SUBROUTINE diffuco_initialize
190
191
192
193	!! ================================================================================================================================
194	!! SUBROUTINE : diffuco_main
195	!!
196	!>\BRIEF The root subroutine for the module, which calls all other required
197	!! subroutines.
198	!!
199	!! DESCRIPTION :
200
201	!! This is the main subroutine for the module.
202	!! First it calculates the surface drag coefficient (via a call to diffuco_aero), using available parameters to determine
203	!! the stability of air in the surface layer by calculating the Richardson Nubmber. If a value for the
204	!! surface drag coefficient is passed down from the atmospheric model and and if the flag 'ldq_cdrag_from_gcm'
205	!! is set to TRUE, then the subroutine diffuco_aerod is called instead. This calculates the aerodynamic coefficient. \n
206	!!
207	!! Following this, an estimation of the saturated humidity at the surface is made (via a call
208	!! to qsatcalc in the module qsat_moisture). Following this the beta coefficients for sublimation (via
209	!! diffuco_snow), interception (diffuco_inter), bare soil (diffuco_bare), and transpiration (via
210	!! diffuco_trans_co2) are calculated in sequence. Finally
211	!! the alpha and beta coefficients are combined (diffuco_comb). \n
212	!!
213	!! The surface drag coefficient is calculated for use within the module enerbil. It is required to to
214	!! calculate the aerodynamic coefficient for every flux. \n
215	!!
216	!! The various beta coefficients are used within the module enerbil for modifying the rate of evaporation,
217	!! as appropriate for the surface. As explained in Chapter 2 of Guimberteau (2010), that module (enerbil)
218	!! calculates the rate of evaporation essentially according to the expression $E = /beta E_{pot}$, where
219	!! E is the total evaporation and $E_{pot}$ the evaporation potential. If $\beta = 1$, there would be
220	!! essentially no resistance to evaporation, whereas should $\beta = 0$, there would be no evaporation and
221	!! the surface layer would be subject to some very stong hydrological stress. \n
222	!!
223	!! The following processes are calculated:
224	!! - call diffuco_aero for aerodynamic transfer coeficient
225	!! - call diffuco_snow for partial beta coefficient: sublimation
226	!! - call diffuco_inter for partial beta coefficient: interception for each type of vegetation
227	!! - call diffuco_bare for partial beta coefficient: bare soil
228	!! - call diffuco_trans_co2 for partial beta coefficient: transpiration for each type of vegetation, using Farqhar's formula
229	!! - call diffuco_comb for alpha and beta coefficient
230	!! - call chemistry_bvoc for alpha and beta coefficients for biogenic emissions
231	!!
232	!! RECENT CHANGE(S): None
233	!!
234	!! MAIN OUTPUT VARIABLE(S): humrel, q_cdrag, vbeta, vbeta1, vbeta4,
235	!! vbeta2, vbeta3, rveget, cimean
236	!!
237	!! REFERENCE(S) :
238	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
239	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
240	!! Circulation Model. Journal of Climate, 6, pp.248-273.
241	!! - de Rosnay, P, 1999. ReprÃ©sentation des interactions sol-plante-atmosphÃšre dans le modÃšle de circulation gÃ©nÃ©rale
242	!! du LMD, 1999. PhD Thesis, UniversitÃ© Paris 6, available (25/01/12):
243	!! http://www.ecmwf.int/staff/patricia_de_rosnay/publications.html#8
244	!! - Ducharne, A, 1997. Le cycle de l'eau: modÃ©lisation de l'hydrologie continentale, Ã©tude de ses interactions avec
245	!! le climat, PhD Thesis, UniversitÃ© Paris 6
246	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
247	!! sur le cycle de l'eau, PhD Thesis, available (25/01/12):
248	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
249	!! - LathiÃšre, J, 2005. Evolution des Ã©missions de composÃ©s organiques et azotÃ©s par la biosphÃšre continentale dans le
250	!! modÃšle LMDz-INCA-ORCHIDEE, UniversitÃ© Paris 6
251	!!
252	!! FLOWCHART :
253	!! \latexonly
254	!! \includegraphics[scale=0.5]{diffuco_main_flowchart.png}
255	!! \endlatexonly
256	!! \n
257	!_ ================================================================================================================================
258
259	SUBROUTINE diffuco_main (kjit, kjpindex, &
260	& index, indexveg, indexlai, u, v, &
261	& zlev, z0m, z0h, roughheight, temp_sol, temp_air, &
262	& temp_growth, rau, q_cdrag, qsurf, qair, &
263	& q2m, t2m, pb, evap_bare_lim, &
264	& evapot, evapot_corr, snow, flood_frac, flood_res, &
265	& frac_nobio, snow_nobio, totfrac_nobio, swnet, swdown, &
266	& coszang, ccanopy, humrel, vegstress, veget, &
267	& veget_max, circ_class_biomass,circ_class_n, qsintveg, qsintmax, assim_param, &
268	& vbeta, vbeta1, vbeta2, vbeta3, &
269	& vbeta3pot, vbeta4, vbeta5, gsmean, rveget, &
270	& rstruct, cimean, gpp, lalo, neighbours, &
271	& resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, &
272	& frac_snow_veg, frac_snow_nobio, &
273	& hist_id, hist2_id, vbeta2sum, vbeta3sum, Isotrop_Abs_Tot_p, &
274	!!$ & Isotrop_Tran_Tot_p, lai_per_level, lai, vbeta23, leaf_ci, &
275	& Isotrop_Tran_Tot_p, lai_per_level, vbeta23, leaf_ci, &
276	& gs_distribution, gs_diffuco_output, gstot_component, gstot_frac,&
277	& warnings, u_speed, profile_vbeta3, profile_rveget, &
278	& delta_c13_assim, leaf_ci_out)
279
280	!! 0. Variable and parameter declaration
281
282	!! 0.1 Input variables
283
284	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
285	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
286	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (-)
287	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier (-)
288	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-)
289	INTEGER(i_std),DIMENSION (kjpindex*(nlai+1)), INTENT (in) :: indexlai !! Indeces of the points on the 3D map
290	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map (-)
291	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
292	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
293	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
294	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0m !! Surface roughness Length for momentum (m)
295	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0h !! Surface roughness Length for heat (m)
296	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective height for roughness (m)
297	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature (K)
298	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature (K)
299	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
300	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
301	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Near surface air specific humidity (kg kg^{-1})
302	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level air specific humidity (kg kg^{-1})
303	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q2m !! 2m air specific humidity (kg kg^{-1})
304	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: t2m !! 2m air temperature (K)
305	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
306	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_frac !! Fraction of floodplains
307	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_res !! Reservoir in floodplains (estimation to avoid over-evaporation)
308	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface level pressure (hPa)
309	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evap_bare_lim !! Limit to the bare soil evaporation when the
310	!! 11-layer hydrology is used (-)
311	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm day^{-1})
312	!! NdN This variable does not seem to be used at
313	!! all in diffuco
314	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation
315	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice,lakes,cities,... (-)
316	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice,lakes,cities,... (kg m^{-2})
317	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+cities+... (-)
318	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2})
319	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux (W m^{-2})
320	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: coszang !! Cosine of the solar zenith angle (unitless)
321	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration inside the canopy (ppm)
322	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (-)
323	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty)
324	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Stand level biomass @tex $(gC m^{-2})$ @endtex
325	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Stand level biomass @tex $(gC m^{-2})$ @endtex
326
327	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2})
328	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception
329	!! (kg m^{-2})
330	REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param !! referenced vcmax, nue, leaf N
331	!! for photosynthesis
332	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
333	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb),INTENT (in):: neighbours !! Vector of neighbours for each
334	!! grid point (1=N, 2=E, 3=S, 4=W)
335	REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: resolution !! The size in m of each grid-box in X and Y
336	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ptnlev1 !! 1st level of soil temperature (Kelvin)
337	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation expressed in mm/tstep
338	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT (in) :: frac_age !! Age efficiency for isoprene emissions (from STOMATE)
339	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
340	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegetated area
341	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in):: frac_snow_nobio !! Snow cover fraction on none vegetated area
342	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: vegstress !! Relative Soil water content accounting for root profile (m3 m-3)
343	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: Isotrop_Abs_Tot_p !! Absorbed radiation per level for photosynthesis
344	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: Isotrop_Tran_Tot_p !! Absorbed radiation per level for photosynthesis
345	REAL(r_std), DIMENSION(:,:,:), INTENT(IN) :: lai_per_level !! This is the LAI per vertical level
346	!! @tex $(m^{2} m^{-2})$
347	REAL(r_std), DIMENSION (:), INTENT (in) :: u_speed !! wind speed at specified canopy level (m/s)
348
349	!! 0.2 Output variables
350
351	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient (-)
352	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation (-)
353	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation (-)
354	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta5 !! Beta for floodplains
355	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gsmean !! Mean stomatal conductance to CO2 (mol m-2 s-1)
356	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss (-)
357	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for transpiration (-)
358	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3pot !! Beta for potential transpiration
359	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Stomatal resistance for the whole canopy (s m^{-1})
360	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: rstruct !! Structural resistance for the vegetation
361	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! Mean leaf Ci (ppm)
362	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gpp !! Assimilation ((gC m^{-2} s^{-1}), total area)
363	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta23 !! Beta for fraction of wetted foliage that will
364	!! transpire once intercepted water has evaporated (-)
365	REAL(r_std),DIMENSION (kjpindex,nvm,nlai), INTENT (out) :: leaf_ci !! intercellular CO2 concentration (ppm)
366	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: gs_distribution
367	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: gs_diffuco_output
368	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: gstot_component
369	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: gstot_frac
370	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: profile_vbeta3
371	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: profile_rveget
372
373	REAL(r_std),DIMENSION (kjpindex,nvm) :: delta_c13_assim !! C13 concentration in delta notation
374	!! @tex $ permille $ @endtex (per thousand)
375	REAL(r_std),DIMENSION (kjpindex,nvm) :: leaf_ci_out !! Ci
376
377	!! 0.3 Modified variables
378
379	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (inout) :: humrel !! Soil moisture stress (within range 0 to 1)
380	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Surface drag coefficient (-)
381	REAL(r_std),DIMENSION(kjpindex,nvm,nwarns), INTENT(inout) &
382	:: warnings !! A warning counter
383
384	!! 0.4 Local variables
385	REAL(r_std),DIMENSION(kjpindex) :: vbeta2sum, vbeta3sum
386	REAL(r_std),DIMENSION(kjpindex) :: raero !! Aerodynamic resistance (s m^{-1})
387	INTEGER(i_std) :: ilai
388	REAL(r_std),DIMENSION(kjpindex,nvm) :: lai !! PFT leaf area index (m^{2} m^{-2})
389	CHARACTER(LEN=4) :: laistring
390	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
391	REAL(r_std),DIMENSION(kjpindex) :: qsatt !! Surface saturated humidity (kg kg^{-1})
392	REAL(r_std),DIMENSION(kjpindex,nvm) :: cim !! Intercellular CO2 over nlai
393	! Isotope_c13
394	INTEGER(i_std), DIMENSION(kjpindex) :: index_assi !! indices (unitless)
395	INTEGER(i_std) :: nia,inia,nina,iainia !! counter/indices (unitless)
396	INTEGER(i_std) :: ji,jl,jv !! Indices (-)
397	REAL(r_std), DIMENSION(nlai+1) :: laitab !! LAI per layer (m^2.m^{-2})
398	REAL(r_std), DIMENSION(kjpindex,nvm) :: assimi !! assimilation (at a specific LAI level)
399	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
400	REAL(r_std), DIMENSION(kjpindex,nvm) :: assimtot !! total assimilation
401	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
402	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot) :: assimi_lev !! assimilation (at all LAI levels)
403	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
404
405	!_ ================================================================================================================================
406
407
408	gs_distribution = 0.0d0
409	wind(:) = SQRT (u(:)u(:) + v(:)v(:))
410
411	!! 1. Calculate the different coefficients
412
413	IF (.NOT.ldq_cdrag_from_gcm) THEN
414	! Case 1a)
415	CALL diffuco_aero (kjpindex, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, &
416	qsurf, qair, snow, q_cdrag)
417	ENDIF
418
419	! Case 1b)
420	CALL diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
421
422	!! 2. Make an estimation of the saturated humidity at the surface
423
424	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
425
426	!! 3. Calculate the beta coefficient for sublimation
427
428	CALL diffuco_snow (kjpindex, qair, qsatt, rau, u, v, q_cdrag, &
429	& snow, frac_nobio, totfrac_nobio, snow_nobio, &
430	& frac_snow_veg, frac_snow_nobio, vbeta1)
431
432
433	CALL diffuco_flood (kjpindex, qair, qsatt, rau, u, v, q_cdrag, evapot, evapot_corr, &
434	& flood_frac, flood_res, vbeta5)
435
436	!! 4. Calculate the beta coefficient for interception
437
438	CALL diffuco_inter (kjpindex, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, &
439	& qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
440
441
442	!! 5. Calculate the beta coefficient for transpiration
443	! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
444	! case 8a)
445	CALL diffuco_trans_co2 (kjpindex, dt_sechiba, swdown, temp_sol, pb, qsurf, q2m, t2m, &
446	temp_growth, rau, u, v, q_cdrag, humrel, vegstress, &
447	assim_param, ccanopy, circ_class_biomass, circ_class_n,&
448	veget_max, qsintveg, qsintmax, vbeta3, vbeta3pot, &
449	rveget, rstruct, cimean, gsmean, gpp, vbeta23, Isotrop_Abs_Tot_p,&
450	Isotrop_Tran_Tot_p, lai_per_level, gs_distribution, gs_diffuco_output, &
451	gstot_component, gstot_frac, veget, warnings, u_speed, profile_vbeta3, &
452	profile_rveget, index_assi, leaf_ci, assimtot, assimi_lev, hist_id, &
453	indexveg, indexlai, index, kjit, cim)
454
455	! calculate C13 isotopic signature
456	IF (ok_c13) THEN
457
458	CALL isotope_c13(kjpindex, index_assi, leaf_ci, ccanopy, lai_per_level, &
459	assimtot, delta_c13_assim, leaf_ci_out, nvm, &
460	nlevels_tot, assimi_lev)
461
462	ENDIF
463
464	!
465	!biogenic emissions
466	!
467	IF ( ok_bvoc ) THEN
468
469	! Calculate BVOC emissions
470	CALL chemistry_bvoc (kjpindex, swdown, coszang, temp_air, &
471	temp_sol, ptnlev1, precip_rain, humrel, veget_max, &
472	circ_class_biomass, circ_class_n, frac_age, lalo, ccanopy, &
473	cim, wind, snow, &
474	veget, hist_id, hist2_id, kjit, index, &
475	indexlai, indexveg)
476
477	ENDIF
478	!
479	! combination of coefficient : alpha and beta coefficient
480	! beta coefficient for bare soil
481	!
482
483	CALL diffuco_bare (kjpindex, evap_bare_lim, vbeta2, vbeta3, vbeta4)
484
485	!! 6. Combine the alpha and beta coefficients
486
487	! Ajout qsintmax dans les arguments de la routine.... Nathalie / le 13-03-2006
488	CALL diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, snow, &
489	veget, circ_class_biomass,circ_class_n, &
490	tot_bare_soil, vbeta1, vbeta2, vbeta3, vbeta4, vbeta, qsintmax, vbeta2sum, vbeta3sum)
491
492	CALL xios_orchidee_send_field("q_cdrag",q_cdrag)
493	CALL xios_orchidee_send_field("raero",raero)
494	CALL xios_orchidee_send_field("wind",wind)
495	CALL xios_orchidee_send_field("qsatt",qsatt)
496	CALL xios_orchidee_send_field("coszang",coszang)
497	CALL xios_orchidee_send_field('cim', cim)
498
499	IF ( .NOT. almaoutput ) THEN
500	CALL histwrite_p(hist_id, 'raero', kjit, raero, kjpindex, index)
501	CALL histwrite_p(hist_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
502	CALL histwrite_p(hist_id, 'Wind', kjit, wind, kjpindex, index)
503	CALL histwrite_p(hist_id, 'qsatt', kjit, qsatt, kjpindex, index)
504
505	IF ( hist2_id > 0 ) THEN
506	CALL histwrite_p(hist2_id, 'raero', kjit, raero, kjpindex, index)
507	CALL histwrite_p(hist2_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
508	CALL histwrite_p(hist2_id, 'Wind', kjit, wind, kjpindex, index)
509	CALL histwrite_p(hist2_id, 'qsatt', kjit, qsatt, kjpindex, index)
510	ENDIF
511	ELSE
512	CALL histwrite_p(hist_id, 'cim', kjit, cim, kjpindex*nvm, indexveg)
513	ENDIF
514
515	IF (printlev>=3) WRITE (numout,*) ' diffuco_main done '
516
517	END SUBROUTINE diffuco_main
518
519	!! =============================================================================================================================
520	!! SUBROUTINE: diffuco_finalize
521	!!
522	!>\BRIEF Write to restart file
523	!!
524	!! DESCRIPTION: This subroutine writes the module variables and variables calculated in diffuco
525	!! to restart file
526	!!
527	!! RECENT CHANGE(S): None
528	!! REFERENCE(S): None
529	!! FLOWCHART: None
530	!! \n
531	!_ ==============================================================================================================================
532	SUBROUTINE diffuco_finalize (kjit, kjpindex, rest_id, rstruct )
533
534	!! 0. Variable and parameter declaration
535	!! 0.1 Input variables
536	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
537	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
538	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (-)
539	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! Structural resistance for the vegetation
540
541	!! 0.4 Local variables
542	INTEGER :: ilai
543	CHARACTER(LEN=4) :: laistring
544	CHARACTER(LEN=80) :: var_name
545
546	!_ ================================================================================================================================
547
548	!! 1. Prepare the restart file for the next simulation
549	IF (printlev>=3) WRITE (numout,*) 'Complete restart file with DIFFUCO variables '
550
551	CALL restput_p (rest_id, 'rstruct', nbp_glo, nvm, 1, kjit, rstruct, 'scatter', nbp_glo, index_g)
552	CALL restput_p (rest_id, 'leaf_ci', nbp_glo, nvm, nlai, kjit, leaf_ci, 'scatter', nbp_glo, index_g)
553
554	END SUBROUTINE diffuco_finalize
555
556
557	!! ================================================================================================================================
558	!! SUBROUTINE : diffuco_clear
559	!!
560	!>\BRIEF Housekeeping module to deallocate the variables
561	!! leaf_ci, rstruct and raero
562	!!
563	!! DESCRIPTION : Housekeeping module to deallocate the variables
564	!! leaf_ci, rstruct and raero
565	!!
566	!! RECENT CHANGE(S) : None
567	!!
568	!! MAIN OUTPUT VARIABLE(S) : None
569	!!
570	!! REFERENCE(S) : None
571	!!
572	!! FLOWCHART : None
573	!! \n
574	!_ ================================================================================================================================
575
576	SUBROUTINE diffuco_clear()
577
578	IF (ALLOCATED (leaf_ci)) DEALLOCATE (leaf_ci)
579
580	! Deallocate and reset variables in chemistry module
581	CALL chemistry_clear
582
583	END SUBROUTINE diffuco_clear
584
585
586	!! ================================================================================================================================
587	!! SUBROUTINE : diffuco_aero
588	!!
589	!>\BRIEF This module first calculates the surface drag
590	!! coefficient, for cases in which the surface drag coefficient is NOT provided by the coupled
591	!! atmospheric model LMDZ or when the flag ldq_cdrag_from_gcm is set to FALSE
592	!!
593	!! DESCRIPTION : Computes the surface drag coefficient, for cases
594	!! in which it is NOT provided by the coupled atmospheric model LMDZ. The module first uses the
595	!! meteorolgical input to calculate the Richardson Number, which is an indicator of atmospheric
596	!! stability in the surface layer. The formulation used to find this surface drag coefficient is
597	!! dependent on the stability determined. \n
598	!!
599	!! Designation of wind speed
600	!! \latexonly
601	!! \input{diffucoaero1.tex}
602	!! \endlatexonly
603	!!
604	!! Calculation of geopotential. This is the definition of Geopotential height (e.g. Jacobson
605	!! eqn.4.47, 2005). (required for calculation of the Richardson Number)
606	!! \latexonly
607	!! \input{diffucoaero2.tex}
608	!! \endlatexonly
609	!!
610	!! \latexonly
611	!! \input{diffucoaero3.tex}
612	!! \endlatexonly
613	!!
614	!! Calculation of the virtual air temperature at the surface (required for calculation
615	!! of the Richardson Number)
616	!! \latexonly
617	!! \input{diffucoaero4.tex}
618	!! \endlatexonly
619	!!
620	!! Calculation of the virtual surface temperature (required for calculation of th
621	!! Richardson Number)
622	!! \latexonly
623	!! \input{diffucoaero5.tex}
624	!! \endlatexonly
625	!!
626	!! Calculation of the squared wind shear (required for calculation of the Richardson
627	!! Number)
628	!! \latexonly
629	!! \input{diffucoaero6.tex}
630	!! \endlatexonly
631	!!
632	!! Calculation of the Richardson Number. The Richardson Number is defined as the ratio
633	!! of potential to kinetic energy, or, in the context of atmospheric science, of the
634	!! generation of energy by wind shear against consumption
635	!! by static stability and is an indicator of flow stability (i.e. for when laminar flow
636	!! becomes turbulent and vise versa). It is approximated using the expression below:
637	!! \latexonly
638	!! \input{diffucoaero7.tex}
639	!! \endlatexonly
640	!!
641	!! The Richardson Number hence calculated is subject to a minimum value:
642	!! \latexonly
643	!! \input{diffucoaero8.tex}
644	!! \endlatexonly
645	!!
646	!! Computing the drag coefficient. We add the add the height of the vegetation to the
647	!! level height to take into account that the level 'seen' by the vegetation is actually
648	!! the top of the vegetation. Then we we can subtract the displacement height.
649	!! \latexonly
650	!! \input{diffucoaero9.tex}
651	!! \endlatexonly
652	!!
653	!! For the stable case (i.e $R_i$ $\geq$ 0)
654	!! \latexonly
655	!! \input{diffucoaero10.tex}
656	!! \endlatexonly
657	!!
658	!! \latexonly
659	!! \input{diffucoaero11.tex}
660	!! \endlatexonly
661	!!
662	!! For the unstable case (i.e. $R_i$ < 0)
663	!! \latexonly
664	!! \input{diffucoaero12.tex}
665	!! \endlatexonly
666	!!
667	!! \latexonly
668	!! \input{diffucoaero13.tex}
669	!! \endlatexonly
670	!!
671	!! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere.
672	!! To prevent this, a minimum limit to the drag coefficient is defined as:
673	!!
674	!! \latexonly
675	!! \input{diffucoaero14.tex}
676	!! \endlatexonly
677	!!
678	!! RECENT CHANGE(S): None
679	!!
680	!! MAIN OUTPUT VARIABLE(S): q_cdrag
681	!!
682	!! REFERENCE(S) :
683	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
684	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
685	!! Circulation Model. Journal of Climate, 6, pp.248-273
686	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
687	!! sur le cycle de l'eau, PhD Thesis, available from:
688	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
689	!! - Jacobson M.Z., Fundamentals of Atmospheric Modeling (2nd Edition), published Cambridge
690	!! University Press, ISBN 0-521-54865-9
691	!!
692	!! FLOWCHART :
693	!! \latexonly
694	!! \includegraphics[scale=0.5]{diffuco_aero_flowchart.png}
695	!! \endlatexonly
696	!! \n
697	!_ ================================================================================================================================
698
699	SUBROUTINE diffuco_aero (kjpindex, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, &
700	qsurf, qair, snow, q_cdrag)
701
702	!! 0. Variable and parameter declaration
703
704	!! 0.1 Input variables
705
706	INTEGER(i_std), INTENT(in) :: kjpindex, kjit !! Domain size
707	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
708	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
709	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first atmospheric layer (m)
710	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0h !! Surface roughness Length for heat (m)
711	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0m !! Surface roughness Length for momentum (m)
712	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective roughness height (m)
713	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Ground temperature (K)
714	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature (K)
715	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! near surface specific air humidity (kg kg^{-1})
716	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
717	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
718
719	!! 0.2 Output variables
720
721	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: q_cdrag !! Surface drag coefficient (-)
722
723	!! 0.3 Modified variables
724
725	!! 0.4 Local variables
726
727	INTEGER(i_std) :: ji, jv
728	REAL(r_std) :: speed, zg, zdphi, ztvd, ztvs, zdu2
729	REAL(r_std) :: zri, cd_neut, zscf, cd_tmp
730	REAL(r_std) :: snowfact
731	!_ ================================================================================================================================
732
733	!! 1. Initialisation
734
735	! test if we have to work with q_cdrag or to calcul it
736	DO ji=1,kjpindex
737
738	!! 1a).1 Designation of wind speed
739	!! \latexonly
740	!! \input{diffucoaero1.tex}
741	!! \endlatexonly
742	speed = wind(ji)
743
744	!! 1a).2 Calculation of geopotentiel
745	!! This is the definition of Geopotential height (e.g. Jacobson eqn.4.47, 2005). (required
746	!! for calculation of the Richardson Number)
747	!! \latexonly
748	!! \input{diffucoaero2.tex}
749	!! \endlatexonly
750	zg = zlev(ji) * cte_grav
751
752	!! \latexonly
753	!! \input{diffucoaero3.tex}
754	!! \endlatexonly
755	zdphi = zg/cp_air
756
757	!! 1a).3 Calculation of the virtual air temperature at the surface
758	!! required for calculation of the Richardson Number
759	!! \latexonly
760	!! \input{diffucoaero4.tex}
761	!! \endlatexonly
762	ztvd = (temp_air(ji) + zdphi / (un + rvtmp2 * qair(ji))) * (un + retv * qair(ji))
763
764	!! 1a).4 Calculation of the virtual surface temperature
765	!! required for calculation of the Richardson Number
766	!! \latexonly
767	!! \input{diffucoaero5.tex}
768	!! \endlatexonly
769	ztvs = temp_sol(ji) * (un + retv * qsurf(ji))
770
771	!! 1a).5 Calculation of the squared wind shear
772	!! required for calculation of the Richardson Number
773	!! \latexonly
774	!! \input{diffucoaero6.tex}
775	!! \endlatexonly
776	zdu2 = MAX(cepdu2,speed**2)
777
778	!! 1a).6 Calculation of the Richardson Number
779	!! The Richardson Number is defined as the ratio of potential to kinetic energy, or, in the
780	!! context of atmospheric science, of the generation of energy by wind shear against consumption
781	!! by static stability and is an indicator of flow stability (i.e. for when laminar flow
782	!! becomes turbulent and vise versa).\n
783	!! It is approximated using the expression below:
784	!! \latexonly
785	!! \input{diffucoaero7.tex}
786	!! \endlatexonly
787	zri = zg * (ztvd - ztvs) / (zdu2 * ztvd)
788
789	!! The Richardson Number hence calculated is subject to a minimum value:
790	!! \latexonly
791	!! \input{diffucoaero8.tex}
792	!! \endlatexonly
793	zri = MAX(MIN(zri,5.),-5.)
794
795	!! 1a).7 Computing the drag coefficient
796	!! We add the add the height of the vegetation to the level height to take into account
797	!! that the level 'seen' by the vegetation is actually the top of the vegetation. Then we
798	!! we can subtract the displacement height.
799	!! \latexonly
800	!! \input{diffucoaero9.tex}
801	!! \endlatexonly
802
803	!! 7.0 Snow smoothering
804	!! Snow induces low levels of turbulence.
805	!! Sensible heat fluxes can therefore be reduced of ~1/3. Pomeroy et al., 1998
806	snowfact=1.
807	IF (snow(ji).GT.snowcri .AND. ok_snowfact .AND. (.NOT. rough_dyn)) snowfact=10.
808
809	cd_neut = ct_karman ** 2. / ( LOG( (zlev(ji) + roughheight(ji)) / z0m(ji) ) * LOG( (zlev(ji) + roughheight(ji)) / z0h(ji) ) )
810
811	!! 1a).7.1 - for the stable case (i.e $R_i$ $\geq$ 0)
812	IF (zri .GE. zero) THEN
813
814	!! \latexonly
815	!! \input{diffucoaero10.tex}
816	!! \endlatexonly
817	zscf = SQRT(un + cd * ABS(zri))
818
819	!! \latexonly
820	!! \input{diffucoaero11.tex}
821	!! \endlatexonly
822	cd_tmp=cd_neut/(un + trois * cb * zri * zscf)
823	ELSE
824
825	!! 1a).7.2 - for the unstable case (i.e. $R_i$ < 0)
826	!! \latexonly
827	!! \input{diffucoaero12.tex}
828	!! \endlatexonly
829	zscf = un / (un + trois * cb * cc * cd_neut * SQRT(ABS(zri) * &
830	& ((zlev(ji) + roughheight(ji)) / z0m(ji)/snowfact)))
831
832	!! \latexonly
833	!! \input{diffucoaero13.tex}
834	!! \endlatexonly
835	cd_tmp=cd_neut * (un - trois * cb * zri * zscf)
836	ENDIF
837
838	!! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere.
839	!! To prevent this, a minimum limit to the drag coefficient is defined as:
840
841	!! \latexonly
842	!! \input{diffucoaero14.tex}
843	!! \endlatexonly
844	!!
845	q_cdrag(ji) = MAX(cd_tmp, 1.e-4/MAX(speed,min_wind))
846
847	! In some situations it might be useful to give an upper limit on the cdrag as well.
848	! The line here should then be uncommented.
849	!q_cdrag(ji) = MIN(q_cdrag(ji), 0.5/MAX(speed,min_wind))
850
851	END DO
852
853	IF (printlev_loc>=3) WRITE (numout,*) ' not ldqcdrag_from_gcm : diffuco_aero done '
854
855	END SUBROUTINE diffuco_aero
856
857
858	!! ================================================================================================================================
859	!! SUBROUTINE : diffuco_snow
860	!!
861	!>\BRIEF This subroutine computes the beta coefficient for snow sublimation.
862	!!
863	!! DESCRIPTION : This routine computes beta coefficient for snow sublimation, which
864	!! integrates the snow on both vegetation and other surface types (e.g. ice, lakes,
865	!! cities etc.) \n
866	!!
867	!! A critical depth of snow (snowcri) is defined to calculate the fraction of each grid-cell
868	!! that is covered with snow (snow/snowcri) while the remaining part is snow-free.
869	!! We also carry out a first calculation of sublimation (subtest) to lower down the beta
870	!! coefficient if necessary (if subtest > snow). This is a predictor-corrector test.
871	!!
872	!! RECENT CHANGE(S): None
873	!!
874	!! MAIN OUTPUT VARIABLE(S): ::vbeta1
875	!!
876	!! REFERENCE(S) :
877	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
878	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
879	!! Circulation Model. Journal of Climate, 6, pp. 248-273
880	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
881	!! sur le cycle de l'eau, PhD Thesis, available from:
882	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
883	!!
884	!! FLOWCHART : None
885	!! \n
886	!_ ================================================================================================================================
887
888	SUBROUTINE diffuco_snow (kjpindex, qair, qsatt, rau, u, v,q_cdrag, &
889	& snow, frac_nobio, totfrac_nobio, snow_nobio, &
890	& frac_snow_veg, frac_snow_nobio, vbeta1)
891
892	!! 0. Variable and parameter declaration
893
894	!! 0.1 Input variables
895
896	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
897	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
898	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1})
899	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3})
900	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
901	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
902	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
903	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg m^{-2})
904	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, cities etc. (-)
905	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice, lakes, cities etc. (-)
906	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice, lakes, cities etc. (-)
907	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
908	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio!! Snow cover fraction on non-vegeted area
909
910
911	!! 0.2 Output variables
912
913	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation (dimensionless ratio)
914
915	!! 0.3 Modified variables
916
917	!! 0.4 Local variables
918
919	REAL(r_std) :: subtest !! Sublimation for test (kg m^{-2})
920	REAL(r_std) :: zrapp !! Modified factor (ratio)
921	REAL(r_std) :: speed !! Wind speed (m s^{-1})
922	REAL(r_std) :: vbeta1_add !! Beta for sublimation (ratio)
923	INTEGER(i_std) :: ji, jv !! Indices (-)
924	!_ ================================================================================================================================
925
926	!! 1. Calculate beta coefficient for snow sublimation on the vegetation\n
927
928	DO ji=1,kjpindex ! Loop over # pixels - domain size
929
930	! Fraction of mesh that can sublimate snow
931	vbeta1(ji) = (un - totfrac_nobio(ji)) * frac_snow_veg(ji)
932
933	! Limitation of sublimation in case of snow amounts smaller
934	! than the atmospheric demand.
935	speed = MAX(min_wind, wind(ji))
936
937	subtest = dt_sechiba * vbeta1(ji) * speed * q_cdrag(ji) * rau(ji) * &
938	& ( qsatt(ji) - qair(ji) )
939
940	IF ( subtest .GT. zero ) THEN
941	zrapp = snow(ji) / subtest
942	IF ( zrapp .LT. un ) THEN
943	vbeta1(ji) = vbeta1(ji) * zrapp
944	ENDIF
945	ENDIF
946
947	END DO ! Loop over # pixels - domain size
948
949	!! 2. Add the beta coefficients calculated from other surfaces types (snow on ice,lakes, cities...)
950
951	DO jv = 1, nnobio ! Loop over # other surface types
952	!!$ !
953	!!$ IF ( jv .EQ. iice ) THEN
954	!!$ !
955	!!$ ! Land ice is of course a particular case
956	!!$ !
957	!!$ DO ji=1,kjpindex
958	!!$ vbeta1(ji) = vbeta1(ji) + frac_nobio(ji,jv)
959	!!$ ENDDO
960	!!$ !
961	!!$ ELSE
962	!
963	DO ji=1,kjpindex ! Loop over # pixels - domain size
964
965	vbeta1_add = frac_nobio(ji,jv) * MAX(MIN(snow_nobio(ji,jv)/snowcri,un), zero)
966
967	! Limitation of sublimation in case of snow amounts smaller than
968	! the atmospheric demand.
969	speed = MAX(min_wind, wind(ji))
970
971	!! Limitation of sublimation by the snow accumulated on the ground
972	!! A first approximation is obtained with the old values of
973	!! qair and qsol_sat: function of temp-sol and pb. (see call of qsatcalc)
974	subtest = dt_sechiba * vbeta1_add * speed * q_cdrag(ji) * rau(ji) * &
975	& ( qsatt(ji) - qair(ji) )
976
977	IF ( subtest .GT. zero ) THEN
978	zrapp = snow_nobio(ji,jv) / subtest
979	IF ( zrapp .LT. un ) THEN
980	vbeta1_add = vbeta1_add * zrapp
981	ENDIF
982	ENDIF
983
984	vbeta1(ji) = vbeta1(ji) + vbeta1_add
985
986	ENDDO ! Loop over # pixels - domain size
987
988	!!$ ENDIF
989
990	ENDDO ! Loop over # other surface types
991
992	IF (printlev>=3) WRITE (numout,*) ' diffuco_snow done '
993
994	END SUBROUTINE diffuco_snow
995
996
997	!! ================================================================================================================================
998	!! SUBROUTINE : diffuco_flood
999	!!
1000	!>\BRIEF This routine computes partial beta coefficient : floodplains
1001	!!
1002	!! DESCRIPTION :
1003	!!
1004	!! RECENT CHANGE(S): None
1005	!!
1006	!! MAIN OUTPUT VARIABLE(S) : vbeta5
1007	!!
1008	!! REFERENCE(S) : None
1009	!!
1010	!! FLOWCHART : None
1011	!! \n
1012	!_ ================================================================================================================================
1013
1014	SUBROUTINE diffuco_flood (kjpindex, qair, qsatt, rau, u, v, q_cdrag, evapot, evapot_corr, &
1015	& flood_frac, flood_res, vbeta5)
1016
1017	! interface description
1018	! input scalar
1019	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1020	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
1021	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity
1022	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
1023	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
1024	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
1025	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
1026	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_res !! water mass in flood reservoir
1027	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_frac !! fraction of floodplains
1028	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Potential evaporation
1029	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr!! Potential evaporation2
1030	! output fields
1031	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta5 !! Beta for floodplains
1032
1033	! local declaration
1034	REAL(r_std) :: subtest, zrapp, speed
1035	INTEGER(i_std) :: ji, jv
1036
1037	!_ ================================================================================================================================
1038	!
1039	! beta coefficient for sublimation for floodplains
1040	!
1041	DO ji=1,kjpindex
1042	!
1043	IF (evapot(ji) .GT. min_sechiba) THEN
1044	vbeta5(ji) = flood_frac(ji) *evapot_corr(ji)/evapot(ji)
1045	ELSE
1046	vbeta5(ji) = flood_frac(ji)
1047	ENDIF
1048	!
1049	! -- Limitation of evaporation in case of water amounts smaller than
1050	! the atmospheric demand.
1051
1052	!
1053	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
1054	!
1055	subtest = dt_sechiba * vbeta5(ji) * speed * q_cdrag(ji) * rau(ji) * &
1056	& ( qsatt(ji) - qair(ji) )
1057	!
1058	IF ( subtest .GT. zero ) THEN
1059	zrapp = flood_res(ji) / subtest
1060	IF ( zrapp .LT. un ) THEN
1061	vbeta5(ji) = vbeta5(ji) * zrapp
1062	ENDIF
1063	ENDIF
1064	!
1065	END DO
1066
1067	IF (printlev>=3) WRITE (numout,*) ' diffuco_flood done '
1068
1069	END SUBROUTINE diffuco_flood
1070
1071
1072	!! ================================================================================================================================
1073	!! SUBROUTINE : diffuco_inter
1074	!!
1075	!>\BRIEF This routine computes the partial beta coefficient
1076	!! for the interception for each type of vegetation
1077	!!
1078	!! DESCRIPTION : We first calculate the dry and wet parts of each PFT (wet part = qsintveg/qsintmax).
1079	!! The former is submitted to transpiration only (vbeta3 coefficient, calculated in
1080	!! diffuco_trans_co2), while the latter is first submitted to interception loss
1081	!! (vbeta2 coefficient) and then to transpiration once all the intercepted water has been evaporated
1082	!! (vbeta23 coefficient). Interception loss is also submitted to a predictor-corrector test,
1083	!! as for snow sublimation. \n
1084	!!
1085	!! \latexonly
1086	!! \input{diffucointer1.tex}
1087	!! \endlatexonly
1088	!! Calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water
1089	!! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss.
1090	!!
1091	!! \latexonly
1092	!! \input{diffucointer2.tex}
1093	!! \endlatexonly
1094	!!
1095	!! Calculation of $\beta_3$, the canopy transpiration resistance
1096	!! \latexonly
1097	!! \input{diffucointer3.tex}
1098	!! \endlatexonly
1099	!!
1100	!! We here determine the limitation of interception loss by the water stored on the leaf.
1101	!! A first approximation of interception loss is obtained using the old values of
1102	!! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc')
1103	!! \latexonly
1104	!! \input{diffucointer4.tex}
1105	!! \endlatexonly
1106	!!
1107	!! \latexonly
1108	!! \input{diffucointer5.tex}
1109	!! \endlatexonly
1110	!!
1111	!! \latexonly
1112	!! \input{diffucointer6.tex}
1113	!! \endlatexonly
1114	!!
1115	!! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction
1116	!! 'zqsvegrap'.
1117	!! \latexonly
1118	!! \input{diffucointer7.tex}
1119	!! \endlatexonly
1120	!!
1121	!! RECENT CHANGE(S): None
1122	!!
1123	!! MAIN OUTPUT VARIABLE(S): ::vbeta2, ::vbeta23
1124	!!
1125	!! REFERENCE(S) :
1126	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1127	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1128	!! Circulation Model. Journal of Climate, 6, pp. 248-273
1129	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1130	!! sur le cycle de l'eau, PhD Thesis, available from:
1131	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1132	!! - Perrier, A, 1975. Etude physique de l'Ã©vaporation dans les conditions naturelles. Annales
1133	!! Agronomiques, 26(1-18): pp. 105-123, pp. 229-243
1134	!!
1135	!! FLOWCHART : None
1136	!! \n
1137	!_ ================================================================================================================================
1138
1139	SUBROUTINE diffuco_inter (kjpindex, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, &
1140	& qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
1141
1142	!! 0 Variable and parameter declaration
1143
1144	!! 0.1 Input variables
1145
1146	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1147	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
1148	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1})
1149	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
1150	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
1151	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
1152	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
1153	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1)
1154	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! vegetation fraction for each type (fraction)
1155	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2})
1156	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2})
1157	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! architectural resistance (s m^{-1})
1158
1159	!! 0.2 Output variables
1160
1161	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss (-)
1162	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta23 !! Beta for fraction of wetted foliage that will
1163	!! transpire (-)
1164
1165	!! 0.4 Local variables
1166
1167	INTEGER(i_std) :: ji, jv !! (-), (-)
1168	REAL(r_std) :: zqsvegrap, ziltest, zrapp, speed !!
1169	!_ ================================================================================================================================
1170
1171	!! 1. Initialize
1172
1173	vbeta2(:,:) = zero
1174	vbeta23(:,:) = zero
1175
1176	!! 2. The beta coefficient for interception by vegetation.
1177
1178	DO jv = 2,nvm
1179
1180	DO ji=1,kjpindex
1181
1182	IF (veget(ji,jv) .GT. min_sechiba .AND. qsintveg(ji,jv) .GT. zero ) THEN
1183
1184	zqsvegrap = zero
1185	IF (qsintmax(ji,jv) .GT. min_sechiba ) THEN
1186
1187	!! \latexonly
1188	!! \input{diffucointer1.tex}
1189	!! \endlatexonly
1190	!!
1191	!! We calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water
1192	!! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss.
1193	!!
1194	zqsvegrap = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv))
1195	END IF
1196
1197	!! \latexonly
1198	!! \input{diffucointer2.tex}
1199	!! \endlatexonly
1200	speed = MAX(min_wind, wind(ji))
1201
1202	!! Calculation of $\beta_3$, the canopy transpiration resistance
1203	!! \latexonly
1204	!! \input{diffucointer3.tex}
1205	!! \endlatexonly
1206	vbeta2(ji,jv) = veget(ji,jv) * zqsvegrap * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv)))
1207
1208	!! We here determine the limitation of interception loss by the water stored on the leaf.
1209	!! A first approximation of interception loss is obtained using the old values of
1210	!! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc')
1211	!! \latexonly
1212	!! \input{diffucointer4.tex}
1213	!! \endlatexonly
1214	ziltest = dt_sechiba * vbeta2(ji,jv) * speed * q_cdrag(ji) * rau(ji) * &
1215	& ( qsatt(ji) - qair(ji) )
1216
1217	IF ( ziltest .GT. zero ) THEN
1218
1219	!! \latexonly
1220	!! \input{diffucointer5.tex}
1221	!! \endlatexonly
1222	zrapp = qsintveg(ji,jv) / ziltest
1223	IF ( zrapp .LT. un ) THEN
1224
1225	!! \latexonly
1226	!! \input{diffucointer6.tex}
1227	!! \endlatexonly
1228	!!
1229	!! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction
1230	!! 'zqsvegrap'.
1231	IF ( humrel(ji,jv) >= min_sechiba ) THEN
1232	vbeta23(ji,jv) = MAX(vbeta2(ji,jv) - vbeta2(ji,jv) * zrapp, zero)
1233	ELSE
1234	! We don't want transpiration when the soil cannot deliver it
1235	vbeta23(ji,jv) = zero
1236	ENDIF
1237
1238	!! \latexonly
1239	!! \input{diffucointer7.tex}
1240	!! \endlatexonly
1241	vbeta2(ji,jv) = vbeta2(ji,jv) * zrapp
1242	ENDIF
1243	ENDIF
1244	END IF
1245	! ! Autre formulation possible pour l'evaporation permettant une transpiration sur tout le feuillage
1246	! !commenter si formulation Nathalie sinon Tristan
1247	! speed = MAX(min_wind, wind(ji))
1248	!
1249	! vbeta23(ji,jv) = MAX(zero, veget(ji,jv) * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv))) - vbeta2(ji,jv))
1250
1251	END DO
1252
1253	END DO
1254
1255	IF (printlev>=3) WRITE (numout,*) ' diffuco_inter done '
1256
1257	END SUBROUTINE diffuco_inter
1258
1259
1260	!! ==============================================================================================================================
1261	!! SUBROUTINE : diffuco_bare
1262	!!
1263	!>\BRIEF This routine computes the partial beta coefficient corresponding to
1264	!! bare soil
1265	!!
1266	!! DESCRIPTION : Bare soil evaporation is submitted to a maximum possible flow (evap_bare_lim) if the 11-layer hydrology is used.\n
1267	!!
1268	!! Calculation of wind speed
1269	!! \latexonly
1270	!! \input{diffucobare1.tex}
1271	!! \endlatexonly
1272	!!
1273	!! The calculation of $\beta_4$
1274	!! \latexonly
1275	!! \input{diffucobare2.tex}
1276	!! \endlatexonly
1277	!!
1278	!! RECENT CHANGE(S): None
1279	!!
1280	!! MAIN OUTPUT VARIABLE(S): ::vbeta4
1281	!!
1282	!! REFERENCE(S) :
1283	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1284	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1285	!! Circulation Model. Journal of Climate, 6, pp.248-273
1286	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1287	!! sur le cycle de l'eau, PhD Thesis, available from:
1288	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1289	!!
1290	!! FLOWCHART : None
1291	!! \n
1292	!_ ================================================================================================================================
1293
1294	SUBROUTINE diffuco_bare (kjpindex, evap_bare_lim, vbeta2, vbeta3, vbeta4)
1295
1296	!! 0. Variable and parameter declaration
1297
1298	!! 0.1 Input variables
1299	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1300	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evap_bare_lim !! limiting factor for bare soil evaporation when the
1301	!! 11-layer hydrology is used (-)
1302	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Beta for Interception
1303	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Beta for Transpiration
1304
1305	!! 0.2 Output variables
1306
1307	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation (-)
1308
1309	!! 0.3 Modified variables
1310
1311	!! 0.4 Local variables
1312	INTEGER(i_std) :: ji
1313	!_ ================================================================================================================================
1314
1315	IF (printlev>=2) WRITE (numout,*) 'Entering diffuco_bare'
1316
1317	!! 1. Calculation of the soil resistance and the beta (beta_4) for bare soil
1318	DO ji = 1, kjpindex
1319
1320	! The limitation by 1-beta2-beta3 is due to the fact that evaporation under vegetation is possible
1321	!! \latexonly
1322	!! \input{diffucobare3.tex}
1323	!! \endlatexonly
1324	vbeta4(ji) = MIN(evap_bare_lim(ji), un - SUM(vbeta2(ji,:)+vbeta3(ji,:)))
1325	END DO
1326
1327	IF (printlev>=3) WRITE (numout,*) ' diffuco_bare done '
1328
1329	END SUBROUTINE diffuco_bare
1330
1331
1332	!! ==============================================================================================================================
1333	!! SUBROUTINE : diffuco_trans_co2
1334	!!
1335	!>\BRIEF This subroutine computes carbon assimilation and stomatal
1336	!! conductance, following respectively Farqhuar et al. (1980) and Ball et al. (1987)
1337	!! and Yin X and Struik PC. 2009.
1338	!!
1339	!! DESCRIPTION :\n
1340	!! *** General:\n
1341	!! The equations are different depending on the photosynthesis mode (C3 versus C4).
1342	!! Based on the work of Farquhar, von Caemmerer and Berry who published in 1980 a
1343	!! biochemical model for C3 photosynthetic rates (the FvCB model), Yin and
1344	!! Struik 2009 propose an analytical algorithm that incorporates a gs model and uses gm as a
1345	!! temperature-dependent parameter for calculating assimilation based on the FvCB
1346	!! model for various environmental scenarios. Yin and Struik 2009 also propose a
1347	!! C4-equivalent version of the FvCB model and an associated analytical algorithm
1348	!! for this model. They believe that applying FvCB-type models to both C3 and C4
1349	!! crops is recommended to accurately predict the response of crop photosynthesis
1350	!! to multiple, interactive environmental variables. \n
1351	!!
1352	!! Besides an analytical solution for photosynthesis the scheme includes a
1353	!! modified Arrhenius function for the temperature dependence that accounts for a
1354	!! decrease of Vcmax and Jmax at high temperatures and a temperature dependent
1355	!! J max/Vcmax ratio (Kattge & Knorr 2007). The temperature response of Vcmax and
1356	!! J max was parametrized with values from reanalysed data in literature (Kattge &
1357	!! Knorr 2007), whereas Vcmax and Jmax at a reference temperature of 25Â° C were
1358	!! derived from observed species-specific values in the TRY database (Kattge et
1359	!! al. 2011).
1360	!!
1361	!! In this version here we replace the old trunk assimilation and conductance
1362	!! which were computed over 20 levels of LAI and then integrated at the canopy level.
1363	!! These artificial LAI levels were introduced in order to allow for a saturation of
1364	!! photosynthesis with LAI, but they have no physical meaning. The approach of
1365	! the two stream radiation transfer model was extended for a multi-layer canopy.
1366	!! The user can choose the numbers of layers in the run.def and the available light at
1367	!! each layers decides about the numbers of layers used for the photosynthesis
1368	!! calculation.
1369	!!
1370	!! NOTE: The variable Isotrop_Abs_Tot_p is the amount of light abosrbed by the canopy,
1371	!! according to the two stream albedo routine. Notice that this is the absorption
1372	!! profile for light coming from an isotropic source (e.g. clouds, and what little
1373	!! light is scattered by the atmosphere on a sunny day). The real absorbed light
1374	!! should have both the incoming radiation from clouds and the incoming radiation
1375	!! from the sun, not just a single number like we have now. If that information
1376	!! ever becomes available, Collim_Abs_Tot should be passed from the albedo routines
1377	!! and used here as well.
1378	!!
1379	!! RECENT CHANGE(S):
1380	!!
1381	!! MAIN OUTPUT VARIABLE(S): beta coefficients, resistances, CO2 intercellular
1382	!! concentration
1383	!!
1384	!! REFERENCE(S) :
1385	!!
1386	!! Yin X, Struik PC. C3 and C4 photosynthesis models: an overview from the
1387	!! perspective of crop modelling. NJAS-Wageningen Journal of Life Sciences 2009.
1388	!! C.J. Bernacchi, A.R. Portis, H. Nakano, S. von Caemmerer, S.P. Long,
1389	!! Temperature response of mesophyll conductance. Implication for the determination
1390	!! of Rubisco enzyme kinetics and for limitations to photosynthesis in vivo, Plant
1391	!! Physiology 130 (2002) 1992â1998.
1392	!! C.J. Bernacchi, E.L. Singsaas, C. Pimentel, A.R. Portis Jr., S.P. Long,
1393	!! Improved temperature response functions for models of Rubisco-limited
1394	!! photosynthesis, Plant, Cell and Environment 24 (2001) 253â259.
1395	!! B.E. Medlyn, E. Dreyer, D. Ellsworth, M. Forstreuter, P.C. Harley, M.U.F.
1396	!! Kirschbaum, X. Le Roux, P. Montpied, J. Strassemeyer, A. Walcroft, K. Wang, D.
1397	!! Loustau, Temperature response of parameters of a biochemically based model of
1398	!! photosynthesis. II. A review of experimental data, Plant, Cell and Environ- ment
1399	!! 25 (2002) 1167â1179.
1400	!! Kattge, J. and Knorr, W.: Temperature acclimation in a biochemical model of
1401	!! photosynthesis: a reanalysis of data from 36 species, Plant Cell Environ., 30
1402	!! (2007) 1176â1190.
1403	!! Keenan, T., et al., Soil water stress and coupled photosynthesisâconductance
1404	!! models: Bridging the gap between conflicting reports on the relative roles of
1405	!! stomatal, mesophyll conductance and biochemical limitations to photosynthesis.
1406	!! Agric. Forest Meteorol. (2010)
1407	!!
1408	!! FLOWCHART : None
1409	!! \n
1410	!_ ================================================================================================================================
1411
1412	SUBROUTINE diffuco_trans_co2 (kjpindex, dt_sechiba, swdown, temp_sol, pb, qsurf, q2m, t2m, &
1413	temp_growth, rau, u, v, q_cdrag, humrel, vegstress, &
1414	assim_param, Ca, circ_class_biomass,circ_class_n,&
1415	veget_max, qsintveg, qsintmax, vbeta3, vbeta3pot, rveget, rstruct, &
1416	cimean, gsmean, gpp, vbeta23, Isotrop_Abs_Tot_p,&
1417	Isotrop_Tran_Tot_p, lai_per_level, gs_distribution, gs_diffuco_output, &
1418	gstot_component, gstot_frac, veget, warnings, u_speed, profile_vbeta3, &
1419	profile_rveget, index_assi, leaf_ci, assimtot, assimi_lev, &
1420	hist_id, indexveg, indexlai, index, kjit, cim)
1421
1422
1423	!
1424	!! 0. Variable and parameter declaration
1425	!
1426
1427	!
1428	!! 0.1 Input variables
1429	!
1430	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1431	REAL(r_std), INTENT (in) :: dt_sechiba !! Time step (s)
1432	REAL(r_std),DIMENSION (:), INTENT (in) :: swdown !! Downwelling short wave flux
1433	!! @tex ($W m^{-2}$) @endtex
1434	REAL(r_std),DIMENSION (:), INTENT (in) :: temp_sol !! Skin temperature (K)
1435	REAL(r_std),DIMENSION (:), INTENT (in) :: pb !! Lowest level pressure (hPa)
1436	REAL(r_std),DIMENSION (:), INTENT (in) :: qsurf !! Near surface specific humidity
1437	!! @tex ($kg kg^{-1}$) @endtex
1438	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q2m !! 2m specific humidity
1439	!! @tex ($kg kg^{-1}$) @endtex
1440	! In off-line mode q2m and qair are the same.
1441	! In off-line mode t2m and temp_air are the same.
1442	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: t2m !! 2m air temperature (K)
1443	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
1444	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! air density @tex ($kg m^{-3}$) @endtex
1445	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
1446	!! @tex ($m s^{-1}$) @endtex
1447	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
1448	!! @tex ($m s^{-1}$) @endtex
1449	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
1450	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (0-1,unitless)
1451	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: assim_param !! min+max+opt temps (K), vcmax, vjmax for
1452	!! photosynthesis
1453	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1454	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: Ca !! CO2 concentration inside the canopy
1455
1456	REAL(r_std),DIMENSION (:,:,:,:,:), INTENT (in) :: circ_class_biomass
1457	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: circ_class_n
1458
1459	!! @tex ($\mu mol mol^{-1}$) @endtex
1460	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vegstress !! Soil moisture stress (0-1,unitless)
1461	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Maximum vegetation fraction of each PFT inside
1462	!! the grid box (0-1, unitless)
1463	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
1464	!! @tex ($kg m^{-2}$) @endtex
1465	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation
1466	!! @tex ($kg m^{-2}$) @endtex
1467	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23 !! Beta for fraction of wetted foliage that will
1468	!! transpire (unitless)
1469	REAL(r_std),DIMENSION (:,:,:),INTENT (in) :: Isotrop_Abs_Tot_p !! Absorbed radiation per level for photosynthesis
1470	REAL(r_std),DIMENSION (:,:,:),INTENT (in) :: Isotrop_Tran_Tot_p !! Transmitted radiation per level for photosynthesis
1471	REAL(r_std),DIMENSION(:,:,:),INTENT (in) :: lai_per_level !! This is the LAI per vertical level
1472	!! @tex $(m^{2} m^{-2})$
1473	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in) :: veget
1474	REAL(r_std), DIMENSION(:), INTENT (in) :: u_speed !! wind speed at specified canopy level (m/s)
1475
1476	INTEGER(i_std),INTENT(in) :: hist_id !! History_ file identifier (-)
1477	INTEGER(i_std),DIMENSION (kjpindex*nvm),INTENT(in) :: indexveg !! Indices of the points on the 3D map (-)
1478	INTEGER(i_std),DIMENSION (kjpindex*(nlevels_tot+1)),INTENT(in) :: indexlai !! Indices of the points on the 3D map
1479	INTEGER(i_std),DIMENSION (kjpindex),INTENT(in) :: index !! Indices of the points on the map (-)
1480	INTEGER(i_std),INTENT(in) :: kjit !! Time step number (-)
1481
1482	!
1483	!! 0.2 Output variables
1484	!
1485	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for Transpiration (unitless)
1486	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3pot !! Beta for Potential Transpiration
1487	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! stomatal resistance of vegetation
1488	!! @tex ($s m^{-1}$) @endtex
1489	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! structural resistance @tex ($s m^{-1}$) @endtex
1490	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! mean intercellular CO2 concentration
1491	!! @tex ($\mu mol mol^{-1}$) @endtex
1492	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gsmean !! mean stomatal conductance to CO2 (umol m-2 s-1)
1493	REAL(r_Std),DIMENSION (kjpindex,nvm), INTENT (out) :: gpp !! Assimilation ((gC m^{-2} s^{-1}), total area)
1494
1495	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1496	:: gs_distribution
1497	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1498	:: gs_diffuco_output
1499	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1500	:: gstot_component
1501	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1502	:: gstot_frac
1503	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1504	:: profile_vbeta3
1505	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1506	:: profile_rveget
1507	INTEGER(i_Std),DIMENSION (kjpindex), INTENT (out) :: index_assi !! pixels for a given PFT where photosynthesis takes place
1508	REAL(r_std),DIMENSION (kjpindex,nvm,nlai), INTENT (out) :: leaf_ci !! intercellular CO2 concentration (ppm)
1509	REAL(r_std),DIMENSION(kjpindex,nvm), INTENT (out) :: assimtot !! total assimilation
1510	REAL(r_std),DIMENSION(kjpindex,nvm,nlevels_tot), INTENT (out) &
1511	:: assimi_lev !! assimilation (at all LAI levels)
1512	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1513
1514	!+++CHECK+++
1515	! INTENT was changed to inout because it is not yet calculated
1516	! should be set back to out once it is calculated (see ORCHIDEE-CN).
1517	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: cim !! Intercellular CO2 over nlai
1518	!+++++++++++
1519
1520	!
1521
1522
1523	!! 0.3 Modified variables
1524	REAL(r_std),DIMENSION(kjpindex,nvm,nwarns), INTENT(inout) &
1525	:: warnings !! A warning counter
1526
1527	!
1528	!! 0.4 Local variables
1529	!
1530
1531	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot+1) :: assimi !! assimilation (at a specific LAI level)
1532	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1533	REAL(r_std), DIMENSION(kjpindex,nvm) :: coupled_lai
1534	REAL(r_std), DIMENSION(kjpindex,nvm) :: total_lai
1535	!
1536	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot) :: gs_store
1537	REAL(r_Std),DIMENSION (kjpindex,nvm) :: gs_column_total
1538	REAL(r_Std),DIMENSION (kjpindex,nvm) :: gstot_store
1539
1540
1541	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_gs_store
1542	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_cimean
1543	REAL(r_std),DIMENSION (kjpindex,nlevels_tot) :: profile_assimtot
1544
1545	REAL(r_std),DIMENSION (kjpindex,nlevels_tot) :: profile_Rdtot
1546	REAL(r_std),DIMENSION (kjpindex,nlevels_tot) :: profile_gstot
1547	REAL(r_std),DIMENSION (kjpindex,nlevels_tot) :: profile_gstop
1548
1549	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_gsmean
1550
1551	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_gpp
1552
1553
1554	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_rstruct
1555	REAL(r_std),DIMENSION (nlevels_tot) :: profile_speed
1556	REAL(r_std),DIMENSION (nlevels_tot) :: profile_cresist
1557
1558	!
1559	!
1560	REAL(r_std),DIMENSION(kjpindex,nvm,nlevels_tot) :: Abs_light !! fraction of absorbed light within each layer (just to check)
1561	REAL(r_std),DIMENSION(kjpindex,nvm) :: Abs_light_can !! fraction of absorbed light within the canopy (just to check)
1562	REAL(r_std),DIMENSION (kjpindex,nvm) :: vcmax !! maximum rate of carboxylation
1563	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
1564	REAL(r_std),DIMENSION (kjpindex,nvm) :: nue !! Nitrogen use Efficiency with impact of leaf age (umol CO2 (gN)-1 s-1)
1565	REAL(r_std),DIMENSION (kjpindex,nvm) :: ntot !! Total canopy (eg leaf) N (gN m-2[ground])
1566	REAL(r_std) :: N_profile !! N-content reduction factor for a specific canopy depth (unitless)
1567	REAL(r_std),DIMENSION (kjpindex) :: leafN_top !! Leaf N content - top of canopy (gN m-2[leaf])
1568	INTEGER(i_std) :: ji, jv, ilev, limit_photo !! Indices (unitless)
1569	REAL(r_std),DIMENSION(kjpindex,nlevels_tot+1) :: info_limitphoto
1570	REAL(r_std),DIMENSION(kjpindex) :: leaf_ci_lowest !! intercellular CO2 concentration at the lowest
1571	!! LAI level
1572	!! @tex ($\mu mol mol^{-1}$) @endtex
1573	REAL(r_std), DIMENSION(kjpindex) :: zqsvegrap !! relative water quantity in the water
1574	!! interception reservoir (0-1,unitless)
1575	REAL(r_std) :: speed !! wind speed @tex ($m s^{-1}$) @endtex
1576	! Assimilation
1577	LOGICAL, DIMENSION(kjpindex) :: assimilate !! where assimilation is to be calculated
1578	!! (unitless)
1579	INTEGER(i_std) :: nia,inia,nina !! counter/indices (unitless)
1580	INTEGER(i_std) :: inina,iainia !! counter/indices (unitless)
1581	INTEGER(i_std), DIMENSION(kjpindex) :: index_non_assi !! indices (unitless)
1582	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: vc2 !! rate of carboxylation (at a specific LAI level)
1583	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
1584
1585	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: vj2 !! rate of Rubisco regeneration (at a specific LAI
1586	!! level) @tex ($\mu mol e- m^{-2} s^{-1}$) @endtex
1587	REAL(r_std), DIMENSION(kjpindex) :: gstop !! stomatal conductance to H2O at topmost level
1588	!! @tex ($m s^{-1}$) @endtex
1589	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: gs !! stomatal conductance to CO2
1590	!! @tex ($\mol m^{-2} s^{-1}$) @endtex
1591	REAL(r_std), DIMENSION(kjpindex,nlevels_tot) :: templeafci
1592
1593
1594	REAL(r_std), DIMENSION(kjpindex) :: gamma_star !! CO2 compensation point (ppm)
1595	!! @tex ($\mu mol mol^{-1}$) @endtex
1596
1597	REAL(r_std), DIMENSION(kjpindex) :: air_relhum !! air relative humidity at 2m
1598	!! @tex ($kg kg^{-1}$) @endtex
1599	REAL(r_std), DIMENSION(kjpindex) :: VPD !! Vapor Pressure Deficit (kPa)
1600	REAL(r_std), DIMENSION(kjpindex) :: water_lim !! water limitation factor (0-1,unitless)
1601
1602	REAL(r_std), DIMENSION(kjpindex) :: gstot !! total stomatal conductance to H2O
1603	!! Final unit is
1604	!! @tex ($m s^{-1}$) @endtex
1605	REAL(r_std), DIMENSION(kjpindex) :: Rdtot !! Total Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 mâ2 sâ1)
1606	REAL(r_std), DIMENSION(kjpindex) :: gs_top !! leaf stomatal conductance to H2O across all top levels
1607	!! @tex ($\mol H2O m^{-2} s^{-1}$) @endtex
1608	REAL(r_std), DIMENSION(kjpindex) :: qsatt !! surface saturated humidity at 2m (??)
1609	!! @tex ($g g^{-1}$) @endtex
1610	!! levels (0-1,unitless)
1611	REAL(r_std), DIMENSION(kjpindex) :: T_Vcmax !! Temperature dependance of Vcmax (unitless)
1612	REAL(r_std), DIMENSION(kjpindex) :: S_Vcmax_acclim_temp !! Entropy term for Vcmax
1613	!! accounting for acclimation to temperature (J K-1 mol-1)
1614	REAL(r_std), DIMENSION(kjpindex) :: T_Jmax !! Temperature dependance of Jmax
1615	REAL(r_std), DIMENSION(kjpindex) :: S_Jmax_acclim_temp !! Entropy term for Jmax
1616	!! accounting for acclimation to temperature (J K-1 mol-1)
1617	REAL(r_std), DIMENSION(kjpindex) :: T_gm !! Temperature dependance of gmw
1618	REAL(r_std), DIMENSION(kjpindex) :: T_Rd !! Temperature dependance of Rd (unitless)
1619	REAL(r_std), DIMENSION(kjpindex) :: T_Kmc !! Temperature dependance of KmC (unitless)
1620	REAL(r_std), DIMENSION(kjpindex) :: T_KmO !! Temperature dependance of KmO (unitless)
1621	REAL(r_std), DIMENSION(kjpindex) :: T_Sco !! Temperature dependance of Sco
1622	REAL(r_std), DIMENSION(kjpindex) :: T_gamma_star !! Temperature dependance of gamma_star (unitless)
1623	REAL(r_std), DIMENSION(kjpindex) :: vc !! Maximum rate of Rubisco activity-limited carboxylation (mumol CO2 mâ2 sâ1)
1624	REAL(r_std), DIMENSION(kjpindex) :: vj !! Maximum rate of e- transport under saturated light (mumol CO2 mâ2 sâ1)
1625	REAL(r_std), DIMENSION(kjpindex) :: gm !! Mesophyll diffusion conductance (molCO2 mÃ¢2 sÃ¢1 barÃ¢1)
1626	REAL(r_std), DIMENSION(kjpindex) :: g0var
1627	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: Rd !! Day respiration (respiratory CO2 release other than by photorespiration)(mumol CO2 mÃ¢2 sÃ¢1)
1628	REAL(r_std), DIMENSION(kjpindex) :: Kmc !! MichaelisâMenten constant of Rubisco for CO2 (mubar)
1629	REAL(r_std), DIMENSION(kjpindex) :: KmO !! MichaelisâMenten constant of Rubisco for O2 (mubar)
1630	REAL(r_std), DIMENSION(kjpindex) :: Sco !! Relative CO2 /O2 specificity factor for Rubisco (bar bar-1)
1631	REAL(r_std), DIMENSION(kjpindex) :: gb_co2 !! Boundary-layer conductance (molCO2 mÃ¢2 sÃ¢1 barÃ¢1)
1632	REAL(r_std), DIMENSION(kjpindex) :: gb_h2o !! Boundary-layer conductance (molH2O mÃ¢2 sÃ¢1 barÃ¢1)
1633	REAL(r_std), DIMENSION(kjpindex) :: fvpd !! Factor for describing the effect of leaf-to-air vapour difference on gs (-)
1634	REAL(r_std), DIMENSION(kjpindex) :: low_gamma_star !! Half of the reciprocal of Sc/o (bar bar-1)
1635	REAL(r_std) :: fcyc !! Fraction of electrons at PSI that follow cyclic transport around PSI (-)
1636	REAL(r_std) :: z !! A lumped parameter (see Yin et al. 2009) ( mol mol-1)
1637	REAL(r_std) :: Rm !! Day respiration in the mesophyll (umol CO2 mâ2 sâ1)
1638	REAL(r_std) :: Cs_star !! Cs -based CO2 compensation point in the absence of Rd (ubar)
1639	REAL(r_std), DIMENSION(kjpindex) :: Iabs !! Photon flux density absorbed by leaf photosynthetic pigments (umol photon mâ2 sâ1)
1640	REAL(r_std), DIMENSION(kjpindex) :: Jmax !! Maximum value of J under saturated light (umol eâ mâ2 sâ1)
1641	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: JJ !! Rate of eâ transport (umol eâ mâ2 sâ1)
1642	REAL(r_std) :: J2 !! Rate of all eâ transport through PSII (umol eâ mâ2 sâ1)
1643	REAL(r_std) :: VpJ2 !! eâ transport-limited PEP carboxylation rate (umol CO2 mâ2 sâ1)
1644	REAL(r_std) :: A_1, A_3 !! Lowest First and third roots of the analytical solution for a general
1645	!! cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 mâ2 sâ1)
1646	REAL(r_std) :: A_1_tmp, A_3_tmp !! Temporary First and third roots of the analytical solution
1647	!! for a general cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 mâ2 sâ1)
1648	REAL(r_std) :: Obs !! Bundle-sheath oxygen partial pressure (ubar)
1649	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: Cc !! Chloroplast CO2 partial pressure (ubar)
1650	REAL(r_std) :: ci_star !! Ci -based CO2 compensation point in the absence of Rd (ubar)
1651	REAL(r_std) :: a,b,c,d,m,f,j,g,h,i,l,p,q,r !! Variables used for solving the cubic equation (see Yin et al. (2009))
1652	REAL(r_std) :: QQ,UU,PSI,x1,x2,x3 !! Variables used for solving the cubic equation (see Yin et al. (2009))
1653	REAL(r_std) :: cresist !! coefficient for resistances (??)
1654	rEAL(r_std) :: lai_min !! mininum lai to calculate photosynthesis
1655
1656	LOGICAL,DIMENSION(kjpindex,nlevels_tot) :: top_level !! A flag to see if a given level is considered
1657	!! a "top" level or not. Must be done this way
1658	!! because there could be several "top" levels.
1659	REAL(r_std) :: totlai_top !! The sum of the LAI in all the levels that we
1660	!! declare as the "top".
1661	INTEGER(i_std) :: lev_count !! How many levels we have already taken for
1662	!! the "top"
1663
1664	! @defgroup Photosynthesis Photosynthesis
1665	! @{
1666	! 1. Preliminary calculations\n
1667	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: light_tran_to_level !! Cumulative amount of light transmitted to a given level
1668	LOGICAL,DIMENSION (kjpindex) :: lskip !! A flag to indicate we've hit an unusual
1669	!! case were we may have numerical instability.
1670
1671	INTEGER :: ic, jc, kc
1672
1673	INTEGER :: ipts
1674	LOGICAL :: llight
1675	LOGICAL :: llai
1676	INTEGER,SAVE :: istep=0, counter=0
1677
1678	!_ ================================================================================================================================
1679
1680	IF (printlev>=4) WRITE (numout,*) ' diffuco_trans_co2'
1681	IF (printlev_loc>=5) WRITE(numout,*) 'diffuco_trans_co2: nlai,lai_per_level',nlai,lai_per_level(:,:,:)
1682
1683	gs_store = 0.0d0
1684	gs_column_total = 0.0d0
1685	gs_distribution = 0.0d0
1686	gs_diffuco_output = 0.0d0
1687
1688	gstot_frac = 0.0d0
1689	gstot_component = 0.d0
1690	gstot_store = 0.0d0
1691
1692	profile_gs_store = 0.0d0
1693	profile_cimean = 0.0d0
1694	profile_assimtot = 0.0d0
1695
1696	profile_Rdtot = 0.0d0
1697	profile_gstot = 0.0d0
1698	profile_gstop = 0.0d0
1699
1700	profile_gsmean = 0.0d0
1701	profile_gpp = 0.0d0
1702	profile_rveget = 2.5d20
1703	profile_rstruct = 0.0d0
1704	profile_speed = 0.0d0
1705	profile_cresist = 0.0d0
1706	profile_vbeta3 = 0.0d0
1707
1708	! Debug
1709	IF(printlev_loc>=4)THEN
1710	WRITE(numout,*) 'Absorbed light in diffuco. ',test_grid,test_pft
1711	DO ilev=nlevels_tot,1,-1
1712	WRITE(numout,'(I5,10F20.10)') &
1713	ilev, Isotrop_Abs_Tot_p(test_grid,test_pft,ilev),&
1714	Isotrop_Tran_Tot_p(test_grid,test_pft,ilev),&
1715	lai_per_level(test_grid,test_pft,ilev)
1716	ENDDO
1717	ENDIF
1718	!-
1719
1720	gs_column_total = 0.0d0
1721	gs_distribution = 0.0d0
1722	gs_diffuco_output = 0.0d0
1723	gstot_frac = 0.0d0
1724	gstot_component = 0.d0
1725	gstot_store = 0.0d0
1726	gs_store(:,:,:)=zero
1727	!
1728	!
1729	! Photosynthesis parameters
1730	!
1731	nue(:,:) = assim_param(:,:,inue)
1732	ntot(:,:) = assim_param(:,:,ileafn)
1733	!vcmax(:,:) = assim_param(:,:,ivcmax) ! from before nitrogen is included
1734	!WRITE(numout,*) '210515 vcmax:', vcmax(:,:)
1735	!
1736	! @addtogroup Photosynthesis
1737	! @{
1738	! 1.3 Estimate relative humidity of air (for calculation of the stomatal conductance).\n
1739	!! \latexonly
1740	!! \input{diffuco_trans_co2_1.3.tex}
1741	!! \endlatexonly
1742	! @}
1743	!
1744	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
1745	air_relhum(:) = &
1746	( qsurf(:) * pb(:) / (Tetens_1+qsurf(:)* Tetens_2) ) / &
1747	( qsatt(:)pb(:) / (Tetens_1+qsatt(:)Tetens_2 ) )
1748
1749	! WRITE(numout, *) '050315 air_relhum is: ', air_relhum
1750
1751	VPD(:) = ( qsatt(:)pb(:) / (Tetens_1+qsatt(:)Tetens_2 ) ) &
1752	- ( qsurf(:) * pb(:) / (Tetens_1+qsurf(:)* Tetens_2) )
1753	! VPD is needed in kPa
1754	! We limit the impact of VPD in the range of [0:6] kPa
1755	! This seems to be the typical range of values
1756	! that the vapor pressure difference can take, although
1757	! there is a paper by Mott and Parkhurst, Plant, Cell & Environment,
1758	! Vol. 14, pp. 509--519 (1991) which gives a value of 15 kPa.
1759	!!$ VPD(:) = MAX(zero,MIN(VPD(:)/10.,6.))
1760	VPD(:) = VPD(:)/10.
1761
1762	!
1763	! 2. beta coefficient for vegetation transpiration
1764	!
1765	rstruct(:,1) = rstruct_const(1)
1766	rveget(:,:) = undef_sechiba
1767	vbeta3(:,:) = zero
1768	vbeta3pot(:,:) = zero
1769	gsmean(:,:) = zero
1770	gpp(:,:) = zero
1771	light_tran_to_level(:,:,:)=un
1772	assimi_lev(:,:,:) = zero
1773	cim(:,:) = zero
1774	cimean(:,1) = Ca(:)
1775
1776	!
1777	! 2.1 Initializations
1778	! Loop over vegetation types
1779	!
1780
1781	DO jv = 2,nvm
1782	gamma_star(:) = zero
1783	Kmo(:) = zero
1784	Kmc(:) = zero
1785	gm(:) = zero
1786	g0var(:) = zero
1787
1788	Cc(:,:) = zero
1789	vc2(:,:) = zero
1790	vj2(:,:) = zero
1791	JJ(:,:) = zero
1792	info_limitphoto(:,:) = zero
1793	gs(:,:) = zero
1794	templeafci(:,:) = zero
1795	assimi(:,:,:) = zero
1796	Rd(:,:) = zero
1797
1798	!
1799	! beta coefficient for vegetation transpiration
1800	!
1801	rstruct(:,jv) = rstruct_const(jv)
1802	cimean(:,jv) = Ca(:)
1803
1804	DO ipts = 1, kjpindex
1805
1806	total_lai(ipts, jv) = cc_to_lai( circ_class_biomass(ipts, jv,:, ileaf,&
1807	icarbon),circ_class_n(ipts, jv,:), jv )
1808
1809	IF( jv .eq. test_pft .and. printlev_loc >= 4 ) &
1810	WRITE(numout,*) '101116 total_lai',total_lai(ipts,jv)
1811
1812	ENDDO
1813
1814	! ---TEMP---
1815	! calculation of vcmax before N is included
1816	! IF (downregulation_co2) THEN
1817	! vcmax(:,jv) = assim_param(:,jv,ivcmax)*&
1818	! (un-downregulation_co2_coeff(jv)*&
1819	! log(Ca(:)/downregulation_co2_baselevel))
1820	! ELSE
1821	! vcmax(:,jv) = assim_param(:,jv,ivcmax)
1822	! ENDIF
1823	!-----------
1824
1825	WHERE (total_lai(:,jv) .GT. min_sechiba)
1826
1827	leafN_top(:) = ntot(:,jv) * ext_coeff_N(jv) / &
1828	( 1. -exp(-ext_coeff_N(jv) * total_lai(:,jv)) )
1829
1830	ELSEWHERE
1831
1832	leafN_top(:) = zero
1833
1834	ENDWHERE
1835
1836	vcmax(:,jv) = nue(:,jv)*leafN_top(:)
1837
1838	!WRITE(numout,*) '210515 assim_param(grid,pft,npco2)','grid:'grid,'pft:',pft,'npco2:',npco2
1839	!WRITE(numout,*) assim_param(:,:,:)
1840	!! mask that contains points where there is photosynthesis
1841	!! For the sake of vectorisation [DISPENSABLE], computations are done only for convenient points.
1842	!! nia is the number of points where the assimilation is calculated and nina the number of points where photosynthesis is not
1843	!! calculated (based on criteria on minimum or maximum values on LAI, vegetation fraction, shortwave incoming radiation,
1844	!! temperature and relative humidity).
1845	!! For the points where assimilation is not calculated, variables are initialized to specific values.
1846	!! The assimilate(kjpindex) array contains the logical value (TRUE/FALSE) relative to this photosynthesis calculation.
1847	!! The index_assi(kjpindex) array indexes the nia points with assimilation, whereas the index_no_assi(kjpindex) array indexes
1848	!! the nina points with no assimilation.
1849	!
1850	! set minimum lai value to test if photosynthesis should be calculated
1851	! 0.0001 is abitrary
1852	! could be externalised
1853	lai_min = 0.0001
1854	!
1855	nia=0
1856	nina=0
1857	index_assi=0
1858	index_non_assi=0
1859	llai=.FALSE.
1860	!
1861	DO ji=1,kjpindex
1862
1863	IF(veget_max(ji,jv) == zero)THEN
1864	! this vegetation type is not present, so no reason to do the
1865	! calculation
1866	CYCLE
1867	ENDIF
1868
1869	!
1870	! This is tricky. We need a certain amount of LAI to have photosynthesis.
1871	! However, we divide our canopy up into certain levels, so even if we
1872	! have a total amount that is sufficient, there might not be enough
1873	! in any particular layer. So we check each layer to see if there
1874	! is one that has enough to give a non-zero photosynthesis. In
1875	! the worst-case scenario (small LAI in every other level), this becomes
1876	! our top level for the gs calculation. If gs_top_level is zero below,
1877	! we have a divide by zero issue.
1878	llai=.FALSE.
1879	DO ilev=1,nlevels_tot
1880	IF(lai_per_level(ji,jv,ilev) .GT. lai_min) THEN
1881	llai=.TRUE.
1882	ELSE
1883	ENDIF
1884	ENDDO
1885
1886
1887	IF ( llai .AND. &
1888	( veget_max(ji,jv) .GT. min_sechiba ) ) THEN
1889
1890	! Here we need to do something different than in the previous model.
1891	! We explicitly need to check for absorbed light. This occurs
1892	! because our albedo in the two stream model is a function
1893	! of the solar angle. If the sun has set, we have no
1894	! light, and therefore we cannot have photosynthesis.
1895	llight=.FALSE.
1896	DO ilev=1,nlevels_tot
1897	IF(Isotrop_Abs_Tot_p(ji,jv,ilev) .GT. min_sechiba) llight=.TRUE.
1898	ENDDO
1899	IF ( ( swdown(ji) .GT. min_sechiba ) .AND. &
1900	( humrel(ji,jv) .GT. min_sechiba) .AND. &
1901	( temp_growth(ji) .GT. tphoto_min(jv) ) .AND. &
1902	llight .AND. &
1903	( temp_growth(ji) .LT. tphoto_max(jv) )) THEN
1904	!
1905	assimilate(ji) = .TRUE.
1906	nia=nia+1
1907	index_assi(nia)=ji
1908	!
1909	ELSE
1910	!
1911	assimilate(ji) = .FALSE.
1912	nina=nina+1
1913	index_non_assi(nina)=ji
1914	!
1915	ENDIF
1916	ELSE
1917	!
1918	assimilate(ji) = .FALSE.
1919	nina=nina+1
1920	index_non_assi(nina)=ji
1921	!
1922	ENDIF
1923	!
1924	ENDDO
1925	!
1926
1927	IF(printlev_loc>=4 .AND. jv == test_pft)THEN
1928	WRITE(numout,*) 'jv: ',jv
1929	WRITE(numout,*) 'test_grid: ',test_grid
1930	WRITE(numout,*) 'llai, veget_max: ',llai,veget_max(:,jv)
1931	WRITE(numout,*) 'temp_growth: ',temp_growth(test_grid)
1932	WRITE(numout,*) 'tphoto_min: ',tphoto_min(jv)
1933	WRITE(numout,*) 'tphoto_max: ',tphoto_max(jv)
1934	WRITE(numout,*) 'swdown: ',swdown(test_grid)
1935	WRITE(numout,*) 'index_non_assi(:): ',index_non_assi(test_grid)
1936	WRITE(numout,*) 'index_assi(:): ',index_assi(test_grid)
1937	WRITE(numout,*) 'nia: ',nia
1938	WRITE(numout,*) 'nina',nina
1939	ENDIF
1940
1941	top_level(:,:)=.FALSE.
1942	gstot(:) = zero
1943	gstop(:) = zero
1944	assimtot(:,jv) = zero
1945	Rdtot(:)=zero
1946	gs_top(:) = zero
1947	lskip(:)=.FALSE.
1948	!
1949	zqsvegrap(:) = zero
1950	WHERE (qsintmax(:,jv) .GT. min_sechiba)
1951	!! relative water quantity in the water interception reservoir
1952	zqsvegrap(:) = MAX(zero, qsintveg(:,jv) / qsintmax(:,jv))
1953	ENDWHERE
1954
1955	!! Calculate the water limitation factor that restricts Vcmax when hydrol_arch is turned off.
1956	! If hydrol_architecture is used, VCmax is not controled but instead there is a control on
1957	! stomatal conductance
1958	IF(ok_hydrol_arch) THEN
1959	WHERE ( assimilate(:) )
1960	water_lim(:) = un
1961	ELSEWHERE
1962	water_lim(:) = zero
1963	ENDWHERE
1964	ELSE
1965	WHERE ( assimilate(:) )
1966	water_lim(:) = MAX( min_sechiba, MIN( vegstress(:,jv)/0.5, un ))
1967	ELSEWHERE
1968	water_lim(:) = zero
1969	ENDWHERE
1970	ENDIF ! (ok_hydrol_arch)
1971
1972	! give a default value of ci for all pixel that do not assimilate
1973	DO ilev=1,nlevels_tot
1974	DO inina=1,nina
1975	leaf_ci(index_non_assi(inina),jv,ilev) = Ca(index_non_assi(inina))
1976	ENDDO
1977	ENDDO
1978
1979	!
1980	!
1981	! Here is the calculation of assimilation and stomatal conductance
1982	! based on the work of Farquahr, von Caemmerer and Berry (FvCB model)
1983	! as described in Yin et al. 2009
1984	! Yin et al. developed a extended version of the FvCB model for C4 plants
1985	! and proposed an analytical solution for both photosynthesis pathways (C3 and C4)
1986	! Photosynthetic parameters used are those reported in Yin et al.
1987	! Except For Vcmax25, relationships between Vcmax25 and Jmax25 for which we use
1988	! Medlyn et al. (2002) and Kattge & Knorr (2007)
1989	! Because these 2 references do not consider mesophyll conductance, we neglect this term
1990	! in the formulations developed by Yin et al.
1991	! Consequently, gm (the mesophyll conductance) tends to the infinite
1992	! This is of importance because as stated by Kattge & Knorr and Medlyn et al.,
1993	! values of Vcmax and Jmax derived with different model parametrizations are not
1994	! directly comparable and the published values of Vcmax and Jmax had to be standardized
1995	! to one consistent formulation and parametrization
1996
1997	! See eq. 6 of Yin et al. (2009)
1998	! Parametrization of Medlyn et al. (2002) - from Bernacchi et al. (2001)
1999	T_KmC(:) = Arrhenius(kjpindex,temp_sol,298.,E_KmC(jv))
2000	T_KmO(:) = Arrhenius(kjpindex,temp_sol,298.,E_KmO(jv))
2001	T_Sco(:) = Arrhenius(kjpindex,t2m,298.,E_Sco(jv))
2002	T_gamma_star(:) = Arrhenius(kjpindex,temp_sol,298.,E_gamma_star(jv))
2003
2004
2005	! Parametrization of Yin et al. (2009) - from Bernacchi et al. (2001)
2006	T_Rd(:) = Arrhenius(kjpindex,temp_sol,298.,E_Rd(jv))
2007
2008
2009	! For C3 plants, we assume that the Entropy term for Vcmax and Jmax
2010	! acclimates to temperature as shown by Kattge & Knorr (2007) - Eq. 9 and 10
2011	! and that Jmax and Vcmax respond to temperature following a modified Arrhenius function
2012	! (with a decrease of these parameters for high temperature) as in Medlyn et al. (2002)
2013	! and Kattge & Knorr (2007).
2014	! In Yin et al. (2009), temperature dependance to Vcmax is based only on a Arrhenius function
2015	! Concerning this apparent unconsistency, have a look to the section 'Limitation of
2016	! Photosynthesis by gm' of Bernacchi (2002) that may provide an explanation
2017
2018	! Growth temperature tested by Kattge & Knorr range from 11 to 35Â°C
2019	! So, we limit the relationship between these lower and upper limits
2020	S_Jmax_acclim_temp(:) = aSJ(jv) + bSJ(jv) * MAX(11., MIN(temp_growth(:),35.))
2021	T_Jmax(:) = Arrhenius_modified(kjpindex,temp_sol,298.,E_Jmax(jv),D_Jmax(jv),S_Jmax_acclim_temp)
2022
2023	S_Vcmax_acclim_temp(:) = aSV(jv) + bSV(jv) * MAX(11., MIN(temp_growth(:),35.))
2024	T_Vcmax(:) = Arrhenius_modified(kjpindex,temp_sol,298.,E_Vcmax(jv),D_Vcmax(jv),S_Vcmax_acclim_temp)
2025
2026	vc(:) = vcmax(:,jv) * T_Vcmax(:)
2027
2028
2029	! As shown by Kattge & Knorr (2007), we make use
2030	! of Jmax25/Vcmax25 ratio (rJV) that acclimates to temperature for C3 plants
2031	! rJV is written as a function of the growth temperature
2032	! rJV = arJV + brJV * T_month
2033	! See eq. 10 of Kattge & Knorr (2007)
2034	! and Table 3 for Values of arJV anf brJV
2035	! Growth temperature is monthly temperature (expressed in Â°C) - See first paragraph of
2036	! section Methods/Data of Kattge & Knorr
2037	vj(:) = ( arJV(jv) + brJV(jv) * MAX(11., MIN(temp_growth(:),35.)) ) * vcmax(:,jv) * T_Jmax(:)
2038
2039	T_gm(:) = Arrhenius_modified(kjpindex,t2m,298.,E_gm(jv),D_gm(jv),S_gm(jv))
2040	gm(:) = gm25(jv) * T_gm(:) * MAX(1-stress_gm(jv), water_lim(:))
2041
2042	g0var(:) = g0(jv)* MAX(1-stress_gs(jv), water_lim(:))
2043	! @endcodeinc
2044	!
2045	KmC(:)=KmC25(jv)*T_KmC(:)
2046	KmO(:)=KmO25(jv)*T_KmO(:)
2047	Sco(:)=Sco25(jv)*T_sco(:)
2048	gamma_star(:) = gamma_star25(jv)*T_gamma_star(:)
2049	! low_gamma_star is defined by Yin et al. (2009)
2050	! as the half of the reciprocal of Sco - See Table 2
2051	!!$ We derive its value from Equation 7 of Yin et al. (2009)
2052	!!$ assuming a constant O2 concentration (Oi)
2053	!!$ low_gamma_star(:) = gamma_star(:)/Oi
2054	low_gamma_star(:) = 0.5 / Sco(:)
2055
2056	! VPD expressed in kPa
2057	! Note : MIN(1.-min_sechiba,MAX(min_sechiba,(a1(jv) - b1(jv) * VPD(:))))
2058	! is always between 0-1 not including 0 and 1
2059	fvpd(:) = 1. / ( 1. / MIN(1.-min_sechiba,MAX(min_sechiba,(a1(jv) - b1(jv) * VPD(:)))) - 1. ) &
2060	* MAX(1-stress_gs(jv), water_lim(:))
2061
2062	! leaf boundary layer conductance
2063	! conversion from a conductance in (m s-1) to (mol H2O m-2 s-1)
2064	! from Pearcy et al. (1991, see below)
2065	gb_h2o(:) = gb_ref * 44.6 * (tp_00/t2m(:)) * (pb(:)/pb_std)
2066
2067	! conversion from (mol H2O m-2 s-1) to (mol CO2 m-2 s-1)
2068	gb_co2(:) = gb_h2o(:) / ratio_H2O_to_CO2
2069
2070	!
2071	! 2.4 Loop over LAI levels to estimate assimilation and conductance
2072	!
2073	IF (nia .GT. 0) THEN
2074
2075	assim_loop:DO inia=1,nia
2076
2077	!
2078	iainia=index_assi(inia)
2079
2080	!
2081	! Find the top layer that still contains lai. We have also
2082	! added a check in here to make sure this layer has absorbed
2083	! light. If there is no absorbed light, we will not calculate
2084	! stomatal conductance for the top level, which means gstop
2085	! will be zero later and we'll crash when we divide by it.
2086	! Include a warning so we know how often this happens.
2087	! The old system of chosing just one top level caused some bad results
2088	! in the transpir, evidently due to the calculation of vbeta3 being
2089	! too low due to too low of LAI in the top level. Therefore, we are
2090	! now going to take several top levels in an effort to make sure
2091	! we have enough LAI. We aim for an arbitrary amount of 0.1 or
2092	! the top three levels. This has not been rigoursly tested, but it
2093	! is theoretically justified by the fact that if the LAI of the top
2094	! level is very low, other leaves will be getting a lot of sun and
2095	! well-exposed to the atmosphere, and therefore contribute equal
2096	! to the resistence of the vegetation.
2097	totlai_top=zero
2098	lev_count=0
2099	!
2100	DO ilev=nlevels_tot, 1, -1
2101	IF(lai_per_level(iainia,jv,ilev) .GT. lai_min)THEN
2102	IF(Isotrop_Abs_Tot_p(iainia,jv,ilev) .GT. min_sechiba)THEN
2103	IF(totlai_top .LT. lai_top(jv) .AND. lev_count .LT. nlev_top)THEN
2104	lev_count=lev_count+1
2105	top_level(iainia,ilev)=.TRUE.
2106	totlai_top=totlai_top+lai_per_level(iainia,jv,ilev)
2107	ELSE
2108	! We have enough levels now.
2109	!WRITE(numout,*) 'Top levels ',totlai_top,lev_count
2110	EXIT
2111	ENDIF
2112	ELSE
2113	IF(err_act.GT.1)THEN
2114	WRITE(numout,*) 'WARNING: Found a level with LAI, ' // &
2115	'but no absorbed light.'
2116	WRITE(numout,*) 'WARNING: iainia,jv,ilev, ',iainia,jv,ilev
2117	WRITE(numout,*) 'WARNING: Isotrop_Abs_Tot_p(iainia,jv,ilev),' //&
2118	'lai_per_level(iainia,jv,ilev), ',iainia,jv,ilev
2119	ENDIF
2120	ENDIF
2121	ENDIF
2122	ENDDO
2123	IF(totlai_top .LT. min_sechiba)THEN
2124	! There seem to be some fringe cases where this happens, usually
2125	! at sunrise or sunset with low LAI levels. Instead of stopping
2126	! here, I'm going to issue a warning. If this case happens,
2127	! we produce no GPP, but we're not going to crash. We need to
2128	! be careful in one of the loops below, since a stomatal
2129	! conductance of zero will cause a crash.
2130	IF(err_act.GT.1)THEN
2131	WRITE(numout,*) 'WARNING: Did not find a level with LAI in diffuco!'
2132	WRITE(numout,*) 'WARNING: iainia,ivm: ',iainia,jv
2133	ENDIF
2134	assimtot(iainia,jv) = zero
2135	Rdtot(iainia) = zero
2136	gstot(iainia) = zero
2137	assimi_lev = zero
2138	lskip(iainia)=.TRUE.
2139	CYCLE assim_loop
2140	ENDIF
2141
2142	! We need to know how much is transmitted to this level
2143	! from the top of the canopy, since this aboslute number
2144	! is what is actually used in the photosyntheis. Right now
2145	! Isotrop_Tran_Tot_p and Isotropic_Abs_Tot_p are the
2146	! relative quantities transmitted through and absorbed
2147	! by the layer.
2148	DO ilev = nlevels_tot-1,1,-1
2149
2150	light_tran_to_level(iainia,jv,ilev) = &
2151	light_tran_to_level(iainia,jv,ilev+1) * &
2152	Isotrop_Tran_Tot_p(iainia,jv,ilev+1)
2153
2154	ENDDO
2155
2156	! Debug
2157	IF(printlev_loc>=4 .AND. jv == test_pft .AND. iainia == test_grid)THEN
2158	WRITE(numout,*) 'Absolute transmitted light to each level',&
2159	jv,iainia
2160	DO ilev = nlevels_tot,1,-1
2161	WRITE(numout,*) ilev,light_tran_to_level(iainia,jv,ilev)
2162	ENDDO
2163	ENDIF
2164	!-
2165
2166	DO ilev = 1, nlevels_tot
2167
2168	IF (Isotrop_Abs_Tot_p(iainia,jv,ilev) .GT. min_sechiba) THEN
2169
2170	IF(lai_per_level(iainia,jv,ilev) .LE. min_sechiba)THEN
2171
2172	! This is not necessarily a problem. LAI is updated in
2173	! slowproc, while absorbed light is updated in condveg,
2174	! which is called after diffuco. So we might be one
2175	! timestep off. It generally should be at night,
2176	! though, so there is no absorbed light. If this
2177	! persists, that is a real problem.
2178	IF(err_act.GT.1)THEN
2179	WRITE(numout,*) 'WARNING: We have absorbed light but no LAI ' // &
2180	'in this level. Skipping.'
2181	WRITE(numout,'(A,3I8)') 'WARNING: iainia,jv,ilev: ',iainia,jv,ilev
2182	WRITE(numout,'(A,2F20.10)') 'WARNING: Isotrop_Abs_Tot_p(iainia,jv,ilev), ' // &
2183	'lai_per_level(iainia,jv,ilev): ',&
2184	Isotrop_Abs_Tot_p(iainia,jv,ilev), &
2185	lai_per_level(iainia,jv,ilev)
2186	ENDIF
2187	CYCLE
2188	ENDIF
2189
2190	! nitrogen is scaled according to the transmitted light
2191	! before the caclulation was: N_profile =
2192	! exp( -ext_coeff_N(jv)*laitab(jl) ). Where laitab is the
2193	! lai contained in a fixed lai layer. ext_coeff_N is
2194	! prescribed so this was a prescribed profile. We replaced
2195	! it by the transmitted light calculated by the 2stream
2196	! model, maybe ( un - Isotrop_Tran_Tot_p(iainia,jv,ilev) )
2197	! could be replaced by the directly calculated absorbed
2198	! light Isotrop_Abs_Tot_p(iainia,jv,ilev).
2199	! Notice that light in the old code was the cumulative
2200	! light transmitted to the layer, while Isotrop_Tran_Tot_p
2201	! is the relative light transmitted through the layer.
2202	!+++TEMP+++
2203	! test for multilayer run
2204	! N_profile = ( un - .7_r_std * &
2205	! ( un - Isotrop_Tran_Tot_p(iainia,jv,ilev) ) )
2206	!++++++++++++++
2207	N_profile = ( un - .7_r_std * &
2208	( un - light_tran_to_level(iainia,jv,ilev) ) )
2209
2210	! When hydrol_arch is turned on there is no water stress on
2211	! the vcmax and photosynthesis is only restricted by
2212	! stomatal closure when transpiration is higher than the
2213	! water supply to the leaves as calculated in hydrol_arch.
2214	! This more consistent with the substantial consensus in
2215	! plant physiology that the main restriction of
2216	! photosynthesis during drought occurs through diffusion
2217	! limitation rather than photochemical or biochemical failure
2218	! (Flexas et al. 2006). Effects of drought stress on
2219	! respiration are less well known: increases, decreases and
2220	! no change have been found. For now there is no effect of
2221	! drought stress on dark respiration when hydrol_arch is
2222	! used.
2223	vc2(iainia,ilev) = vc(iainia) * N_profile * MAX(1-stress_vcmax(jv), water_lim(iainia))
2224	vj2(iainia,ilev) = vj(iainia) * N_profile * MAX(1-stress_vcmax(jv), water_lim(iainia))
2225
2226	! see Comment in legend of Fig. 6 of Yin et al. (2009)
2227	! Rd25 is assumed to equal 0.01 Vcmax25
2228	Rd(iainia,ilev) = vcmax(iainia,jv) * N_profile * 0.01 * &
2229	T_Rd(iainia) * MAX(1-stress_vcmax(jv), water_lim(iainia))
2230	Iabs(iainia) =swdown(iainia) * W_to_mol * RG_to_PAR *&
2231	Isotrop_Abs_Tot_p(iainia,jv,ilev) * &
2232	light_tran_to_level(iainia,jv,ilev) / &
2233	lai_per_level(iainia,jv,ilev)
2234
2235	! If we have a very dense canopy, the lower levels might
2236	! get almost no light. If the light is close to zero, it
2237	! causes a numerical problem below, as x1 will be zero. Of
2238	! course, if the absorbed light is really small
2239	! photosynthesis will not happen, so we can just skip this
2240	! level. This number is fairly arbitrary, but it should be
2241	! as small as possible.
2242	IF(Iabs(iainia) .LT. 1e-12)CYCLE
2243
2244	! eq. 4 of Yin et al (2009)
2245	! Note that Iabs is very sensitive to the number of layers
2246	! use in the photosynthesis module and to the canopy
2247	! structure. The first sensitivity is related to
2248	! none-linearities to the light conditions and is an
2249	! numerical issue. The sensitivity to the canopy structure
2250	! is a wanted behavior (that's why canopy structure was
2251	! implemented) but makes a substantial difference for the
2252	! total GPP in tropical compared to temperate forests. An
2253	! easy way to prescribe a different canopy structure is by
2254	! changing the range of the prescribed height. The parameter
2255	! ::alpha_LL describes the relationship between Jmax and
2256	! light saturation. It is thought to depend on the amount of
2257	! chlorophyl in the leaves which in turn depends on the
2258	! N-content. For the moment we use a prescribed value. After
2259	! the introduction of the structured canopy we noticed an
2260	! important sensitivity of GPP to this parameter which may
2261	! indicate that a substantial part of our canopy is under
2262	! light saturation.
2263	Jmax(iainia) = vj2(iainia,ilev)
2264	JJ(iainia,ilev) = ( alpha_LL(jv) * Iabs(iainia) + Jmax(iainia) - &
2265	sqrt((alpha_LL(jv) * Iabs(iainia) + Jmax(iainia) )**2. &
2266	- 4 * theta(jv) * Jmax(iainia) * alpha_LL(jv) * &
2267	Iabs(iainia)) ) &
2268	/ ( 2 * theta(jv))
2269
2270
2271	IF (printlev_loc>=4 .AND. test_pft == jv) THEN
2272	WRITE(numout,) '************************'
2273	WRITE(numout,*) 'jv,ilev: ',jv,ilev
2274	WRITE(numout,*) 'JJ(iainia,ilev),Jmax(iainia),Iabs(iainia),Rd(iainia,ilev): ',&
2275	JJ(iainia,ilev),Jmax(iainia),Iabs(iainia),Rd(iainia,ilev)
2276	WRITE(numout,*) 'vj2(iainia,ilev),vc2(iainia,ilev), N_profile:',&
2277	vj2(iainia,ilev),vc2(iainia,ilev),N_profile
2278	WRITE(numout,*) 'Isotrop_Abs_Tot_p(iainia,jv,ilev),Isotrop_Tran_Tot_p(iainia,jv,ilev):',&
2279	Isotrop_Abs_Tot_p(iainia,jv,ilev), &
2280	Isotrop_Tran_Tot_p(iainia,jv,ilev)
2281
2282	WRITE(numout,*) 'water_lim(iainia),T_Rd(iainia):',&
2283	water_lim(iainia),T_Rd(iainia)
2284	WRITE(numout,*) 'vc(iainia),vj(iainia), vcmax(iainia,jv):',&
2285	vc(iainia),vj(iainia),vcmax(iainia,jv)
2286	WRITE(numout,) '************************'
2287	ENDIF
2288
2289	IF ( is_c4(jv) ) THEN
2290	!
2291	! @addtogroup Photosynthesis
2292	! @{
2293	!
2294	! 2.4.2 Assimilation for C4 plants (Collatz et al., 1992)\n
2295	!! \latexonly
2296	!! \input{diffuco_trans_co2_2.4.2.tex}
2297	!! \endlatexonly
2298	! @}
2299	!
2300	! Analytical resolution of the Assimilation based
2301	! Yin et al. (2009)
2302	! Eq. 28 of Yin et al. (2009)
2303	fcyc= 1. - ( 4.(1.-fpsir(jv))(1.+fQ(jv)) + &
2304	3.h_protons(jv)fpseudo(jv) ) / &
2305	( 3.h_protons(jv) - 4.(1.-fpsir(jv)))
2306
2307	! See paragraph after eq. (20b) of Yin et al.
2308	Rm=Rd(iainia,ilev)/2.
2309
2310	! We assume that cs_star equals ci_star
2311	! (see Comment in legend of Fig. 6 of Yin et al. (2009)
2312	! Equation 26 of Yin et al. (2009)
2313	Cs_star = (gbs(jv) * low_gamma_star(iainia) * &
2314	Oi - ( 1. + low_gamma_star(iainia) * alpha(jv) / &
2315	0.047) * Rd(iainia,ilev) + Rm ) / ( gbs(jv) + kp(jv) )
2316
2317	! eq. 11 of Yin et al (2009)
2318	J2 = JJ(iainia,ilev) / ( 1. - fpseudo(jv) / ( 1. - fcyc ) )
2319
2320	! Equation right after eq. (20d) of Yin et al. (2009)
2321	z = ( 2. + fQ(jv) - fcyc ) / ( h_protons(jv) * (1. - fcyc ))
2322
2323	VpJ2 = fpsir(jv) * J2 * z / 2.
2324
2325	A_3=9999.
2326
2327	! See eq. right after eq. 18 of Yin et al. (2009)
2328	! The loop over limit_photo calculates the assimultion
2329	! rate of the Rubisco-limited CO2 assimilation (Ac) and
2330	! the rate of electron transport-limited CO2 assimultion
2331	! (Aj). The overall assimilation will be the slower
2332	! of these two processes.
2333	DO limit_photo=1,2
2334	! Is Vc limiting the Assimilation
2335	IF ( limit_photo .EQ. 1 ) THEN
2336	a = 1. + kp(jv) / gbs(jv)
2337	b = 0.
2338	x1 = vc2(iainia,ilev)
2339	x2 = KmC(iainia)/KmO(iainia)
2340	x3 = KmC(iainia)
2341	! Is J limiting the Assimilation
2342	ELSE
2343	a = 1.
2344	b = VpJ2
2345	x1 = (1.- fpsir(jv)) * J2 * z / 3.
2346	x2 = 7. * low_gamma_star(iainia) / 3.
2347	x3 = 0.
2348	ENDIF
2349
2350	m=fvpd(iainia)-g0var(iainia)/gb_co2(iainia)
2351	d=g0var(iainia)(Ca(iainia)-Cs_star) + fvpd(iainia)Rd(iainia,ilev)
2352	f=(b-Rm-low_gamma_star(iainia)Oigbs(jv))x1d + agbs(jv)x1Ca(iainia)d
2353	j=(b-Rm+gbs(jv)x3 + x2gbs(jv)Oi)m + (alpha(jv)x2/0.047-1.)d &
2354	+ agbs(jv)(Ca(iainia)*m - d/gb_co2(iainia) - (Ca(iainia) - Cs_star ))
2355
2356	g=(b-Rm-low_gamma_star(iainia)Oigbs(jv))x1m - (alpha(jv)low_gamma_star(iainia)/0.047+1.)x1*d &
2357	+ agbs(jv)x1(Ca(iainia)m - d/gb_co2(iainia) - (Ca(iainia)-Cs_star ))
2358
2359	h=-((alpha(jv)low_gamma_star(iainia)/0.047+1.)x1m + (agbs(jv)x1(m-1.))/gb_co2(iainia) )
2360	i= ( b-Rm + gbs(jv)x3 + x2gbs(jv)Oi )d + agbs(jv)Ca(iainia)*d
2361	l= ( alpha(jv)x2/0.047 - 1.)m - (agbs(jv)(m-1.))/gb_co2(iainia)
2362
2363	p = (j-(h-l*Rd(iainia,ilev))) / l
2364	q = (i+j*Rd(iainia,ilev)-g) / l
2365	r = -(f-i*Rd(iainia,ilev)) / l
2366
2367	! See Yin et al. (2009) and Baldocchi (1994)
2368	QQ = ( (p*2._r_std) - 3._r_std q) / 9._r_std
2369	UU = ( 2._r_std* (p*3._r_std) - 9._r_std pq + 27._r_std r) /54._r_std
2370
2371
2372	IF (QQ .GE. 0._r_std) THEN
2373	IF (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) THEN
2374	PSI = ACOS(UU/(QQ**1.5_r_std))
2375	A_3_tmp = -2._r_std * SQRT(QQ) * COS(( PSI + 4._r_std * PI)/3._r_std ) - p / 3._r_std
2376	IF (( A_3_tmp .LT. A_3 )) THEN
2377	A_3 = A_3_tmp
2378	info_limitphoto(iainia,ilev)=2.
2379	ELSE
2380	! In case, J is not limiting the assimilation
2381	! we have to re-initialise a, b, x1, x2 and x3 values
2382	! in agreement with a Vc-limited assimilation
2383	a = 1. + kp(jv) / gbs(jv)
2384	b = 0.
2385	x1 = vc2(iainia,ilev)
2386	x2 = KmC(iainia)/KmO(iainia)
2387	x3 = KmC(iainia)
2388	info_limitphoto(iainia,ilev)=1.
2389	ENDIF
2390	ENDIF
2391	ELSE
2392	WRITE(numout,*) 'WARNING: unexpected value for QQ (diffuco.f90)'
2393	ENDIF
2394
2395	IF ( ( A_3 .EQ. 9999. ) .OR. ( A_3 .LT. (-Rd(iainia,ilev)) ) ) THEN
2396	WRITE(numout,*) 'We have a problem in diffuco_trans_co2'
2397	WRITE(numout,*) 'no real positive solution found for pft:',jv
2398	WRITE(numout,*) 'temp_sol:',temp_sol(iainia)
2399	WRITE(numout,*) 'vpd:',vpd(iainia)
2400	WRITE(numout,*) 'Rd:',Rd(iainia,ilev)
2401	A_3 = -Rd(iainia,ilev)
2402	ENDIF
2403
2404	assimi(iainia,jv,ilev) = A_3
2405
2406	IF(printlev_loc>=4 .AND. jv == test_pft .AND. iainia == test_pft)THEN
2407	WRITE(numout,*) '( x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) ', &
2408	( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) ))
2409	WRITE(numout,*) 'x1,assimi(iainia,jv,ilev),Rd(iainia,ilev): ', &
2410	x1,assimi(iainia,jv,ilev),Rd(iainia,ilev)
2411	WRITE(numout,*) 'jv,ilev,iainia: ', jv,ilev,iainia
2412	ENDIF
2413
2414	!++++ TEMP +++++
2415	IF(( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) )) == zero)THEN
2416	! This causes a problem, since we really should not
2417	! be dividing by zero. However, it seems to happen
2418	! fairly rarely, so all we'll do is cycle to the
2419	! next loop and keep track of this warning so we
2420	! can look at it afterwards.
2421	warnings(iainia,jv,iwphoto)=warnings(iainia,jv,iwphoto)+un
2422
2423	IF(err_act.GT.1)THEN
2424	WRITE(numout,*) 'Diffuco, Getting ready to divide by zero! C4'
2425	WRITE(numout,*) '( x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) ', &
2426	( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) ))
2427	WRITE(numout,*) 'x1,assimi(iainia,jv,ilev),Rd(iainia,ilev): ', &
2428	x1,assimi(iainia,jv,ilev),Rd(iainia,ilev)
2429	WRITE(numout,*) 'jv,ilev,iainia: ', jv,ilev,iainia
2430	ENDIF
2431
2432	CYCLE
2433
2434	ENDIF
2435
2436	! Eq. 24 of Yin et al. (2009)
2437	Obs = ( alpha(jv) * assimi(iainia,jv,ilev) ) / ( 0.047 * gbs(jv) ) + Oi
2438	! Eq. 23 of Yin et al. (2009)
2439	Cc(iainia,ilev) = ( ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) * &
2440	( x2 * Obs + x3 ) + low_gamma_star(iainia) * Obs * x1 ) &
2441	/ MAX(min_sechiba, x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ))
2442	! Eq. 22 of Yin et al. (2009)
2443	leaf_ci(iainia,jv,ilev) = ( Cc(iainia,ilev) - ( b - assimi(iainia,jv,ilev) - Rm ) / gbs(jv) ) / a
2444
2445	! Eq. 25 of Yin et al. (2009)
2446	! It should be Cs instead of Ca but it seems that
2447	! other equations in Appendix C make use of Ca
2448	IF ( ABS( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) .LT. min_sechiba ) THEN
2449	gs(iainia,ilev) = g0var(iainia)
2450	ELSE
2451	gs(iainia,ilev) = g0var(iainia) + ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) / &
2452	( Ca(iainia) - Cs_star ) * fvpd(iainia)
2453	ENDIF
2454	!
2455	ENDDO !! ENDDO about limit photo
2456	!
2457	ELSE
2458	!
2459	! @addtogroup Photosynthesis
2460	! @{
2461	!
2462	! 2.4.3 Assimilation for C3 plants (Farqhuar et al., 1980)\n
2463	!! \latexonly
2464	!! \input{diffuco_trans_co2_2.4.3.tex}
2465	!! \endlatexonly
2466	! @}
2467	!
2468	!
2469	!
2470	A_1=9999.
2471
2472	! See eq. right after eq. 18 of Yin et al. (2009)
2473	DO limit_photo=1,2
2474	! Is Vc limiting the Assimilation
2475	IF ( limit_photo .EQ. 1 ) THEN
2476
2477	IF (printlev_loc>=4) THEN
2478	WRITE(numout,*) 'IF - A_1, Do in limit_photo limit_photo .EG. 1 '
2479	ENDIF
2480	x1 = vc2(iainia,ilev)
2481	! It should be O not Oi (comment from Vuichard)
2482	x2 = KmC(iainia) * ( 1. + 2gamma_star(iainia)Sco(iainia) / KmO(iainia) )
2483	! Is J limiting the Assimilation
2484	ELSE
2485	x1 = JJ(iainia,ilev)/4.
2486	x2 = 2. * gamma_star(iainia)
2487	ENDIF
2488
2489
2490	! See Appendix B of Yin et al. (2009)
2491	a = g0var(iainia) * ( x2 + gamma_star(iainia) ) + &
2492	( g0var(iainia) / gm(iainia) + fvpd(iainia) ) * ( x1 - Rd(iainia,ilev) )
2493	b = Ca(iainia) * ( x1 - Rd(iainia,ilev) ) - gamma_star(iainia) * x1 - Rd(iainia,ilev) * x2
2494	c = Ca(iainia) + x2 + ( 1./gm(iainia) + 1./gb_co2(iainia) ) * ( x1 - Rd(iainia,ilev) )
2495	d = x2 + gamma_star(iainia) + ( x1 - Rd(iainia,ilev) ) / gm(iainia)
2496	m = 1./gm(iainia) + ( g0var(iainia)/gm(iainia) + fvpd(iainia) ) * ( 1./gm(iainia) + 1./gb_co2(iainia) )
2497
2498	p = - ( d + (x1 - Rd(iainia,ilev) ) / gm(iainia) + a * (1./gm(iainia) + 1./gb_co2(iainia) ) + &
2499	( g0var(iainia)/gm(iainia) + fvpd(iainia)) * c) / m
2500
2501	q = ( d * ( x1 - Rd(iainia,ilev) ) + ac + ( g0var(iainia)/gm(iainia) + fvpd(iainia) ) b ) / m
2502	r = - a * b / m
2503
2504	! See Yin et al. (2009)
2505	QQ = ( (p*2._r_std) - 3._r_std q) / 9._r_std
2506	UU = ( 2._r_std* (p*3._r_std) - 9._r_std pq + 27._r_std r) /54._r_std
2507
2508	IF (QQ .GE. 0._r_std) THEN
2509
2510	IF (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) THEN
2511	PSI = ACOS(UU/(QQ**1.5_r_std))
2512	A_1_tmp = -2._r_std * SQRT(QQ) * COS( PSI / 3._r_std ) - p / 3._r_std
2513
2514	IF (( A_1_tmp .LT. A_1 )) THEN
2515	A_1 = A_1_tmp
2516	info_limitphoto(iainia,ilev)=2.
2517	ELSE
2518	! In case, J is not limiting the assimilation
2519	! we have to re-initialise x1 and x2 values
2520	! in agreement with a Vc-limited assimilation
2521	x1 = vc2(iainia,ilev)
2522	! It should be O not Oi (comment from Vuichard)
2523	x2 = KmC(iainia) * ( 1. + 2gamma_star(iainia)Sco(iainia) / KmO(iainia) )
2524
2525	ENDIF
2526	ENDIF
2527	ENDIF
2528	ENDDO
2529
2530	IF ( (A_1 .EQ. 9999.) .OR. ( A_1 .LT. (-Rd(iainia,ilev)) ) ) THEN
2531	WRITE(numout,*) 'We have a problem in diffuco_trans_co2'
2532	WRITE(numout,*) 'no real positive solution found for pft:',jv
2533	WRITE(numout,*) 'temp_sol:',temp_sol(iainia)
2534	WRITE(numout,*) 'vpd:',vpd(iainia)
2535	WRITE(numout,*) 'Setting the solution to -Rd(iainia,ilev):',-Rd(iainia,ilev)
2536	A_1 = -Rd(iainia,ilev)
2537	ENDIF
2538	assimi(iainia,jv,ilev) = A_1
2539
2540
2541	! Eq. 18 of Yin et al. (2009)
2542
2543	IF(printlev_loc>=4 .AND. jv == test_pft .AND. iainia == test_pft)THEN
2544	WRITE(numout,*) '( x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) ', &
2545	( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) ))
2546	WRITE(numout,*) 'x1,assimi(iainia,jv,ilev),Rd(iainia,ilev): ', &
2547	x1,assimi(iainia,jv,ilev),Rd(iainia,ilev)
2548	WRITE(numout,*) 'jv,ilev,iainia: ', jv,ilev,iainia
2549	ENDIF
2550
2551	!++++ TEMP +++++
2552	IF(( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) )) == zero)THEN
2553
2554	! This causes a problem, since we really should not
2555	! be dividing by zero. However, it seems to happen
2556	! fairly rarely, so all we'll do is cycle to the
2557	! next loop and keep track of this warning so we
2558	! can look at it afterwards.
2559	warnings(iainia,jv,iwphoto)=warnings(iainia,jv,iwphoto)+un
2560
2561	IF(err_act.GT.1)THEN
2562	WRITE(numout,*) 'Diffuco, Getting ready to divide by zero! C3'
2563	WRITE(numout,*) '( x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) ', &
2564	( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) ))
2565	WRITE(numout,*) 'x1,assimi(iainia,jv,ilev),Rd(iainia,ilev): ', &
2566	x1,assimi(iainia,jv,ilev),Rd(iainia,ilev)
2567	WRITE(numout,*) 'jv,ilev,iainia: ', jv,ilev,iainia
2568	ENDIF
2569
2570	CYCLE
2571
2572	ENDIF
2573	!++++++++++++++++++++++++++++++
2574
2575	Cc(iainia,ilev) = ( gamma_star(iainia) * x1 + ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) * x2 ) &
2576	/ MAX( min_sechiba, x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) )
2577	! Eq. 17 of Yin et al. (2009)
2578	leaf_ci(iainia,jv,ilev) = Cc(iainia,ilev) + assimi(iainia,jv,ilev) / gm(iainia)
2579	! See eq. right after eq. 15 of Yin et al. (2009)
2580	ci_star = gamma_star(iainia) - Rd(iainia,ilev) / gm(iainia)
2581	!
2582	! Eq. 15 of Yin et al. (2009)
2583	IF ( ABS( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) .LT. min_sechiba ) THEN
2584	gs(iainia,ilev) = g0var(iainia)
2585	ELSE
2586	gs(iainia,ilev) = g0var(iainia) + ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) / ( leaf_ci(iainia,jv,ilev) &
2587	- ci_star ) * fvpd(iainia)
2588	ENDIF
2589
2590	ENDIF !C3 vs C4
2591	!
2592	!
2593	! this is only for output to netcdf
2594	gs_diffuco_output(iainia,jv,ilev) = gs(iainia,ilev)
2595	!
2596	!
2597	! @addtogroup Photosynthesis
2598	! @{
2599	!
2600	!! 2.4.4 Estimatation of the stomatal conductance (Ball et al., 1987).\n
2601	!! \latexonly
2602	!! \input{diffuco_trans_co2_2.4.4.tex}
2603	!! \endlatexonly
2604	! @}
2605	!
2606	!
2607	! keep stomatal conductance of the levels which we decide
2608	! are the top levels. Notice that we multiply by the LAI
2609	! in the level here, so this is no longer a per leaf quantity,
2610	! but the whole quantity.
2611	IF ( top_level(iainia,ilev) ) THEN
2612	gs_top(iainia) = gs_top(iainia)+gs(iainia,ilev)*lai_per_level(iainia,jv,ilev)
2613	!
2614	ENDIF
2615	!
2616	! @addtogroup Photosynthesis
2617	! @{
2618	!
2619	!! 2.4.5 Integration at the canopy level\n
2620	!! \latexonly
2621	!! \input{diffuco_trans_co2_2.4.5.tex}
2622	!! \endlatexonly
2623	! @}
2624	! total assimilation and conductance
2625	!
2626	!++++++++++++++++++++
2627	! The equations give the amount of carbon produced per m**2 if
2628	! leaf surface. Therefore, since we want GPP in units of m**2 of land surface,
2629	! we have to multiply by the LAI.
2630
2631	assimtot(iainia,jv) = assimtot(iainia,jv) + &
2632	assimi(iainia,jv,ilev) * lai_per_level(iainia,jv,ilev)
2633
2634	IF (.NOT. ok_impose_can_structure ) THEN
2635	profile_assimtot(iainia,ilev) = profile_assimtot(iainia,ilev) + &
2636	& assimi(iainia,jv,ilev) * lai_per_level(iainia,jv,ilev)
2637	profile_Rdtot(iainia,ilev) = profile_Rdtot(iainia,ilev) + &
2638	& Rd(iainia,ilev) * lai_per_level(iainia,jv,ilev)
2639	profile_gstot(iainia,ilev) = profile_gstot(iainia,ilev) + &
2640	& gs(iainia,ilev) * lai_per_level(iainia,jv,ilev)
2641	END IF ! (ok_impose_can_structure ) THEN
2642
2643	! assimi(iainia)
2644
2645	Rdtot(iainia) = Rdtot(iainia) + &
2646	Rd(iainia,ilev) * lai_per_level(iainia,jv,ilev)
2647	gstot(iainia) = gstot(iainia) + &
2648	gs(iainia,ilev) * lai_per_level(iainia,jv,ilev)
2649	assimi_lev(iainia,jv,ilev) = assimi(iainia,jv,ilev)
2650
2651	IF(test_pft == jv .AND. printlev_loc>=4)THEN
2652	WRITE(numout,'(A,2I5,10ES20.8)') 'assimi middle',ilev,&
2653	iainia,assimi(iainia,jv,ilev),assimtot(iainia,jv),gs(iainia,ilev),&
2654	gstot(iainia),lai_per_level(iainia,jv,ilev)
2655	ENDIF
2656
2657	!! (JR 100714) we need to retain the gs (stomatal conductance) per level and per PFT for use
2658	!! outside of this routine, in the case that the supply term constraint is applied and the
2659	!! stomatal conductance profile needs to be re-calculated.
2660	gs_store(iainia, jv, ilev) = gs(iainia,ilev)
2661	!
2662
2663	IF (.NOT. ok_impose_can_structure ) THEN
2664	profile_gstot(iainia, ilev) = profile_gstot(iainia, ilev) + &
2665	& gs(iainia,ilev) * lai_per_level(iainia,jv,ilev)
2666	END IF ! (ok_impose_can_structure ) THEN
2667
2668	!! (JR 100714) we also retain the gstot (stomatal conductance over all canopy levels, constrained
2669	!! by LAI, again for the recalculation of the stomatal conductance profile.
2670	gstot_store(iainia, jv) = gstot_store(iainia, jv) + (gs(iainia,ilev) * lai_per_level(iainia,jv,ilev))
2671
2672	!
2673	ENDIF ! if there is absorbed light
2674
2675	ENDDO ! loop over levels
2676
2677	ENDDO assim_loop ! loop over assimulation points
2678
2679	IF(jv==testpft) THEN
2680	templeafci(:,:)=leaf_ci(:,testpft,:)
2681	CALL histwrite_p(hist_id, 'Cc', kjit, Cc, kjpindex*(nlevels_tot+1), indexlai)
2682	CALL histwrite_p(hist_id, 'Vc', kjit, Vc2, kjpindex*(nlevels_tot+1), indexlai)
2683	CALL histwrite_p(hist_id, 'Vj', kjit, JJ, kjpindex*(nlevels_tot+1), indexlai)
2684	CALL histwrite_p(hist_id, 'limitphoto', kjit, info_limitphoto, kjpindex*(nlevels_tot+1), indexlai)
2685	CALL histwrite_p(hist_id, 'gammastar', kjit, gamma_star, kjpindex, index)
2686	CALL histwrite_p(hist_id, 'Kmo', kjit, Kmo, kjpindex, index)
2687	CALL histwrite_p(hist_id, 'Kmc', kjit, Kmc, kjpindex, index)
2688	CALL histwrite_p(hist_id, 'gm', kjit, gm, kjpindex, index)
2689	CALL histwrite_p(hist_id, 'gs', kjit, gs, kjpindex*(nlevels_tot+1), indexlai)
2690	! CALL histwrite_p(hist_id, 'leafci', kjit, templeafci, kjpindex*(nlevels_tot), indexlai)
2691	! CALL histwrite_p(hist_id, 'assimi', kjit, assimi, kjpindexnvm(nlevels_tot+1), indexlai)
2692	CALL histwrite_p(hist_id, 'Rd', kjit, Rd, kjpindex*(nlevels_tot+1), indexlai)
2693	ENDIF
2694
2695	!
2696	!! 2.5 Calculate resistances
2697
2698	DO inia=1,nia
2699
2700	!
2701	iainia=index_assi(inia)
2702	!
2703
2704	! We hit a fringe case above where we could not
2705	! calculate the GPP. Skip this loop so we don't
2706	! divide by zero.
2707	IF(lskip(iainia)) CYCLE
2708
2709
2710	!! Mean stomatal conductance for CO2 (mol m-2 s-1)
2711	gsmean(iainia,jv) = gstot(iainia)
2712	!
2713	! cimean is the "mean ci" calculated in such a way that assimilation
2714	! calculated in enerbil is equivalent to assimtot
2715	!
2716
2717	IF (.NOT. ok_impose_can_structure ) THEN
2718	DO ilev = 1, nlevels_tot
2719	profile_gsmean(iainia, jv, ilev) = profile_gstot(iainia, ilev)
2720	IF(printlev_loc >= 4) WRITE(numout,*) '110215 profile_gsmean by level: ', &
2721	ilev, ' is: ', profile_gsmean(iainia, jv, ilev)
2722	ENDDO ! ilev = 1, nlevels_tot
2723	END IF ! (ok_impose_can_structure ) THEN
2724
2725	IF((assimtot(iainia,jv)+Rdtot(iainia)) .GT. min_sechiba)THEN
2726	cimean(iainia,jv) = (gsmean(iainia,jv)-g0var(iainia)) &
2727	/ (fvpd(iainia)*(assimtot(iainia,jv)+Rdtot(iainia))) &
2728	+ gamma_star(iainia)
2729	ELSE
2730	! According to Nicolas, this variable is completely diagnostic and not
2731	! used for anything. If the analytical photosynthesis is unable to find
2732	! a solution for the A_1 coefficients in any of the levels, assimtot will
2733	! be equal to -Rd, which means we will be dividing by zero. This loop
2734	! is to prevent that, and since cimean is diagnostic is doesn't matter
2735	! what the value is.
2736	cimean(iainia,jv) = val_exp
2737	ENDIF
2738
2739
2740
2741	IF (.NOT. ok_impose_can_structure ) THEN
2742	DO ilev = 1, nlevels_tot
2743	IF((profile_assimtot(iainia,ilev)+profile_Rdtot(iainia,ilev)) .GT. min_sechiba)THEN
2744	profile_cimean(iainia,jv,ilev) = (profile_gs_store(iainia,jv,ilev)-g0(jv)) &
2745	/ (fvpd(iainia)*(profile_assimtot(iainia,ilev)+profile_Rdtot(iainia,ilev))) &
2746	+ gamma_star(iainia)
2747	IF(printlev_loc >= 4) THEN
2748	WRITE(numout,*) '210515 fvpd:', fvpd(iainia)
2749	WRITE(numout,*) 'profile_assimtot:', profile_assimtot(iainia,ilev)
2750	WRITE(numout,*) 'profile Rdtot:', profile_Rdtot(iainia,ilev)
2751	ENDIF ! IF(printlev >= 4)
2752	ELSE
2753	! According to Nicolas, this variable is completely diagnostic and not
2754	! used for anything. If the analytical photosynthesis is unable to find
2755	! a solution for the A_1 coefficients in any of the levels, assimtot will
2756	! be equal to -Rd, which means we will be dividing by zero. This loop
2757	! is to prevent that, and since cimean is diagnostic is doesn't matter
2758	! what the value is.
2759	profile_cimean(iainia,jv,ilev) = val_exp
2760	ENDIF
2761	! WRITE(numout, *) '100215b profile_cimean by level: ', ilev, ' is: ', profile_cimean(iainia, jv, ilev)
2762	ENDDO ! ilev = 1, nlevels_tot
2763	END IF ! (ok_impose_can_structure ) THEN
2764
2765
2766	! conversion from umol m-2 (PFT) s-1 to gC m-2 (mesh area) tstep-1
2767	gpp(iainia,jv) = assimtot(iainia,jv)12e-6veget_max(iainia,jv)*dt_sechiba
2768	istep=istep+1
2769
2770	IF (printlev_loc>=4) THEN
2771	WRITE(numout,*) 'GPP DIFFUCO: ',istep,gpp(iainia,jv),iainia,jv,assimtot(iainia,jv),veget_max(iainia,jv),dt_sechiba
2772	ENDIF
2773
2774	IF (.NOT. ok_impose_can_structure ) THEN
2775	DO ilev = 1, nlevels_tot
2776	profile_gpp(iainia,jv,ilev) = profile_assimtot(iainia, ilev)12e-6veget_max(iainia,jv)*dt_sechiba
2777	IF(printlev_loc >= 4) WRITE(numout, *) '100215b profile_gpp by level: ', &
2778	ilev, ' is: ', profile_gpp(iainia, jv, ilev)
2779	END DO ! ilev = 1, nlevels_tot
2780	END IF ! (ok_impose_can_structure ) THEN
2781
2782	! As in Pearcy, Schulze and Zimmermann
2783	! Measurement of transpiration and leaf conductance
2784	! Chapter 8 of Plant Physiological Ecology
2785	! Field methods and instrumentation, 1991
2786	! Editors:
2787	!
2788	! Robert W. Pearcy,
2789	! James R. Ehleringer,
2790	! Harold A. Mooney,
2791	! Philip W. Rundel
2792	!
2793	! ISBN: 978-0-412-40730-7 (Print) 978-94-010-9013-1 (Online)
2794	!
2795	!
2796	!
2797	!++++++++ CHECK +++++++++
2798	! Do we need to normalise gs_top here by the LAI in the level?
2799	! In the trunk, it seems like they use an LAI of 0.1, so renormalise
2800	! to that. The default value of lai_top is 0.1.
2801	!
2802	!! (JR 100714)
2803	!! gs_top is normalised in the trunk over the entire canopy LAI, as it
2804	!! meant to provide backwards compatability with the 'old' method for calculating
2805	!! the stomatal conductance (in diffuco_trans) that calculates gs according to
2806	!! incident SW radiation (rather than CO2 concentration, as here in diffuco_trans_co2)
2807	!!
2808	!! rstruct in diffuco_trans is an extra resistance term added to reflect the fact
2809	!! that less light reaches the lower canopy, and the total resistance used to calculate
2810	!! vbeta3, and hence the transpiration, is (rstruct + rveget)
2811	!!
2812	!! The rstruct term makes the code more difficult to follow, and is no longer has
2813	!! a physical requirement. We are getting small values of gstop as calculated below,
2814	!! which was resulting in big resistances when rstruct = min_sechiba and
2815	!! rveget = un/gstop and total resistance = (rstruct + rveget) - see commented out
2816	!! code below.
2817	!!
2818	!! The most transparent structure is to set total resistance = rveget = (1 / gstot)
2819	!! This is analogous to the way in which assimtot, and hence canopy GPP, is calculated
2820	!! above
2821
2822	gs_top(iainia)=gs_top(iainia)/totlai_top*lai_top(jv)
2823	!+++++++++++++++++++++++++
2824
2825	gstot(iainia) = mol_to_m_1 (temp_sol(iainia)/tp_00)&
2826	(pb_std/pb(iainia))gstot(iainia)ratio_H2O_to_CO2
2827	gstop(iainia) = mol_to_m_1 * (temp_sol(iainia)/tp_00)*&
2828	(pb_std/pb(iainia))gs_top(iainia)ratio_H2O_to_CO2
2829
2830
2831	IF (.NOT. ok_impose_can_structure ) THEN
2832	DO ilev = 1, nlevels_tot
2833
2834	profile_gsmean(iainia,jv,ilev) = profile_gsmean(iainia,jv,ilev) * 1e-6
2835
2836	profile_gstot(iainia,ilev) = mol_to_m_1 (temp_sol(iainia)/tp_00)&
2837	(pb_std/pb(iainia))profile_gstot(iainia,ilev)ratio_H2O_to_CO2
2838
2839	profile_gstop(iainia,ilev) = mol_to_m_1 * (temp_sol(iainia)/tp_00)*&
2840	(pb_std/pb(iainia))gs_top(iainia)ratio_H2O_to_CO2*&
2841	lai_per_level(iainia,jv,ilev)
2842
2843	END DO ! ilev = 1, nlevels_tot
2844	END IF ! (ok_impose_can_structure ) THEN
2845
2846
2847	!IF (jv == test_pft .AND. test_grid == 1) THEN
2848	! WRITE(numout, '(A20)' ) 'Without water stress:'
2849	! WRITE(numout, '(A20, I5)' ) 'pft type is: ', jc
2850	! WRITE(numout, '(A20, ES15.5)') 'Old_Gs_total:', gstot(iainia)
2851	!ENDIF
2852
2853	! THIS PART OF CORRECTION OF ::gstot SHOULD BE CONSISTENCY WITH
2854	! "SECHIBA_HYDROL_ARCH" where the ::gstot is recalculated, again.
2855
2856	!++++++ CHECK ++++++
2857	! We introduce ::coupled_lai to calculate the :: vbeta3
2858	! In an closed canopy not all the leaves fully interact with the
2859	! atmosphere because leaves can shelter each other. In a more open
2860	! canopy most leaves can interact with the atmosphere. We use the
2861	! ratio of lai over ::Pgap as a proxy for the amount of canopy that
2862	! can interact with the atmosphere. Clearly that amount of lai that
2863	! interacts should never exceed the:: total_lai.
2864	coupled_lai(iainia, jv) = biomass_to_coupled_lai( SUM(circ_class_biomass(iainia, &
2865	jv,:, ileaf, icarbon)*circ_class_n(iainia,jv,:)), veget(iainia,jv), jv)
2866
2867	total_lai(iainia, jv) = cc_to_lai( circ_class_biomass(iainia, jv,:, ileaf,&
2868	icarbon),circ_class_n(iainia, jv,:), jv )
2869
2870	!---TEMP---
2871	!!$ IF (printlev_loc>=4) THEN
2872	!!$ WRITE(numout,*) 'coupled LAI, ', jv, coupled_lai(iainia, jv)
2873	!!$ WRITE(numout,*) 'total LAI, ', jv, total_lai(iainia, jv)
2874	!!$ WRITE(numout,*) 'veget, ', jv, veget(iainia, jv)
2875	!!$ WRITE(numout,*) 'biomass(ileaf), ',biomass(iainia, jv, ileaf, icarbon)
2876	!!$ ENDIF
2877	!----------
2878
2879	!+++++CHECK++++
2880	!Here, we try to use coupled_lai to reset/refine the total stomata conductance
2881	!which interacted with the atmosphere for the transpiration for each pft.
2882	gstot(iainia) = gstot(iainia)* (coupled_lai(iainia,jv)/total_lai(iainia,jv))
2883	!++++++++++++++
2884
2885	!+++++++++CHECK++++++++
2886	!To get the resonable transpiration we had to set :: lai_top to
2887	!::lai
2888	!so we were using the entire canopy. Therefore we no longer use
2889	!::lai_top and simply make use of the entire canopy to calculate the
2890	!vegetation resistance from inverse of ::gstot so, we calculate the
2891	!:: rveget using :: 1/gstot
2892	rveget(iainia, jv) = un / gstot(iainia)
2893	!+++++++++++++++++++++++++++++++++++++++
2894
2895	!rveget(iainia,jv) = un/gstop(iainia)
2896	!!$ rveget_min(iainia,jv) = (defc_plus / kzero(jv)) * &
2897	!!$ (un / SUM(lai_per_level(iainia,jv,:)))
2898	!
2899	! rstruct is the difference between rtot (=1./gstot) and rveget
2900	!
2901	! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif
2902	!rstruct(iainia,jv) = un/gstot(iainia) - &
2903	! rveget(iainia,jv)
2904
2905	!++++++++ CHECK +++++++
2906	!It seems that the calcualtion of rstruct is duplicated with the calculation of rvegt
2907	!So, we set the ::rstruct as zero
2908	!rstruct(iainia,jv) = MAX( un/gstot(iainia) - &
2909	! rveget(iainia,jv), min_sechiba)
2910	!
2911	!
2912	!! wind is a global variable of the diffuco module.
2913	speed = MAX(min_wind, wind(iainia))
2914	!
2915	! beta for transpiration
2916	!
2917	! Corrections Nathalie - 28 March 2006 - on advices of Fred Hourdin
2918	!! Introduction of a potentiometer rveg_pft to settle the rveg+rstruct sum problem in the coupled mode.
2919	!! rveg_pft=1 in the offline mode. rveg_pft is a global variable declared in the diffuco module.
2920	!vbeta3(iainia,jv) = veget_max(iainia,jv) * &
2921	! (un - zqsvegrap(iainia)) * &
2922	! (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + &
2923	! rstruct(iainia,jv))))
2924	!! Global resistance of the canopy to evaporation
2925	!cresist=(un / (un + speed * q_cdrag(iainia) * &
2926	! veget(iainia,jv)/veget_max(iainia,jv) * &
2927	! (rveg_pft(jv)*(rveget(iainia,jv) + rstruct(iainia,jv)))))
2928
2929	! (JR 100714)
2930	! previously the resistance term in this expression was (rveget + rstruct), but we now
2931	! use rveget alone
2932
2933	cresist=(un / (un + speed * q_cdrag(iainia) * &
2934	(rveg_pft(jv)*(rveget(iainia,jv)))))
2935
2936	!+++CHECK+++
2937	!For the interception we assume -as earlier- that the surface area
2938	!that intercepts is well estimated by ::veget. This is wrong because
2939	!every leav can contribute to the interception. Depending on weather
2940	!the parameters for interception are bulk values or single leave. Our
2941	!approach is correct or needd to be improved.
2942	!
2943	vbeta3(iainia,jv) = veget_max(iainia,jv) * &
2944	(un - zqsvegrap(iainia)) * cresist + &
2945	MIN( vbeta23(iainia,jv), veget(iainia,jv) * &
2946	zqsvegrap(iainia) * cresist )
2947	!+++++++++++++++++++++++++
2948
2949	!vbeta3(iainia,jv) = veget_max(iainia,jv) * &
2950	! (un - zqsvegrap(iainia)) * cresist + &
2951	! MIN( vbeta23(iainia,jv), veget_max(iainia,jv) * &
2952	! zqsvegrap(iainia) * cresist )
2953
2954
2955
2956	IF (.NOT. ok_impose_can_structure ) THEN
2957	DO ilev = 1, nlevels_tot
2958
2959	! Fix the gstot equal zero condition to aviod the issue of divided by zero
2960	IF ( profile_gstot(iainia,ilev) .LE. min_sechiba ) THEN
2961	profile_gstot(iainia,ilev) = 1.0d-6
2962	ENDIF
2963
2964	profile_rveget(iainia,jv,ilev) = un/profile_gstot(iainia,ilev)
2965
2966	! need to resolve the rstruct issue
2967	profile_rstruct(iainia,jv,ilev) = MAX (un/profile_gstot(iainia,ilev) - &
2968	& profile_rveget(iainia,jv,ilev), min_sechiba)
2969
2970	profile_speed(ilev) = MAX(min_wind, u_speed(ilev)) ! n.b. need to replace this with another wind profile parameterisation
2971
2972	profile_cresist(ilev) = (un / (un + u_speed(ilev) * q_cdrag(iainia) * &
2973	& (rveg_pft(jv) * (profile_rveget(iainia,jv,ilev) + profile_rstruct(iainia,jv,ilev)) ) ) )
2974
2975	profile_vbeta3(iainia,jv,ilev) = veget_max(iainia,jv) * (un - zqsvegrap(iainia)) * profile_cresist(ilev) + &
2976	& MIN( vbeta23(iainia,jv), &
2977	veget_max(iainia,jv) * zqsvegrap(iainia) * profile_cresist(ilev) )
2978
2979	END DO ! ilev = 1, nlevels_tot
2980	END IF ! (ok_impose_can_structure ) THEN
2981
2982
2983	IF (printlev_loc>=4) THEN
2984	DO ilev = 1, nlevels_tot
2985	WRITE(numout, *) 'profile_rveget(iainia, jv, ilev) is: ', profile_rveget(iainia, jv, ilev)
2986	END DO
2987
2988	DO ilev = 1, nlevels_tot
2989	WRITE(numout, *) 'profile_vbeta3(iainia, jv, ilev) is: ', profile_vbeta3(iainia, jv, ilev)
2990	END DO
2991
2992	DO ilev = 1, nlevels_tot
2993	WRITE(numout, *) 'profile_cresist(ilev) is: ', (profile_cresist(ilev))
2994	END DO
2995	END IF !(printlev_loc)
2996
2997	! vbeta3pot for computation of potential transpiration (needed for irrigation)
2998	vbeta3pot(iainia,jv) = MAX(zero, veget_max(iainia,jv) * cresist)
2999	!
3000	!
3001
3002	ENDDO !loop over grid points
3003
3004	!
3005	ENDIF ! if nia LT zero
3006
3007	ENDDO ! loop over vegetation types
3008
3009
3010	DO ic = 1, kjpindex
3011	DO jc = 1, nvm
3012	!++++++TEMP++++++++++++++++++++++++++++++++++++
3013	! IF (printlev_loc>=4 .AND. jc == test_pft) THEN
3014	! WRITE(714, '(A20, ES20.4)') 'Rcanopy,', (un / (un + speed q_cdrag(ic) &
3015	! (rveg_pft(jc)*(rveget(ic,jc) + rstruct(ic,jc)))))
3016
3017	! WRITE(numout,'(A20, ES20.4)') 'Rcanopy',(un / (un + speed q_cdrag(ic) &
3018	! (rveg_pft(jc)*(rveget(ic,jc) + rstruct(ic,jc)))))
3019	! ENDIF
3020	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3021
3022
3023
3024	IF(veget_max(ic,jc) == zero)THEN
3025	! this vegetation type is not present, so no reason to do the
3026	! calculation
3027	CYCLE
3028	ENDIF
3029
3030	DO kc = 1, nlevels_tot
3031	gs_column_total(ic,jc) = gs_column_total(ic,jc) + gs_store(ic,jc,kc) * lai_per_level(ic,jc,kc)
3032
3033	gstot_component(ic,jc,kc) = gs_store(ic,jc,kc) * lai_per_level(ic,jc,kc)
3034
3035	IF (gstot_store(ic,jc) .eq. 0.0d0) THEN
3036	gstot_frac(ic,jc,kc) = 0.0d0
3037	ELSE
3038	gstot_frac(ic,jc,kc) = gstot_component(ic,jc,kc) / gstot_store(ic,jc)
3039	END IF
3040	END DO ! kc = 1, nlevels_tot
3041	END DO ! jc = 1, nvm
3042	END DO ! ic = 1, kjpindex
3043
3044
3045	Abs_light(:,:,:) = zero
3046	Abs_light_can(:,:) = zero
3047
3048	DO ic = 1, kjpindex
3049	DO jc = 1, nvm
3050	IF(veget_max(ic,jc) == zero)THEN
3051	! this vegetation type is not present, so no reason to do the
3052	! calculation
3053	CYCLE
3054	ENDIF
3055
3056	DO kc = 1, nlevels_tot
3057	IF (gs_column_total(ic,jc) .ne. zero) THEN
3058	gs_distribution(ic,jc,kc) = gs_store(ic,jc,kc)* lai_per_level(ic,jc,kc) / gs_column_total(ic,jc)
3059	ELSE
3060	gs_distribution(ic,jc,kc) = zero
3061	END IF
3062	Abs_light(ic,jc,kc)=Isotrop_Abs_Tot_p(ic,jc,kc)*light_tran_to_level(ic,jc,kc)
3063	END DO ! kc = 1, nlevels_tot
3064	IF (jv==test_pft .AND. printlev_loc>=4) THEN
3065	Abs_light_can(:,:)=SUM(Abs_light(:,:,:),3)
3066	WRITE(numout,*) 'Absorbed light in the canopy', Abs_light_can(:,jc)
3067	ENDIF
3068	END DO ! jc = 1, nvm
3069	END DO ! ic = 1, kjpindex
3070
3071	!!$ !++++++ WRITING OUT SOME STUFF
3072	!!$ counter=counter+1
3073	!!$ WRITE(123,'(I8,20F20.10)') counter,gs_distribution(test_grid,test_pft,:)
3074	!!$ WRITE(124,'(I8,20F20.10)') counter,gstot_frac(test_grid,test_pft,:)
3075	!!$ WRITE(125,'(I8,20F20.10)') counter,gstot_store(test_grid,test_pft)
3076	!!$ WRITE(126,'(I8,20F20.10)') counter,gs_store(test_grid,test_pft,:)
3077	!!$ WRITE(127,'(I8,20F20.10)') counter,lai_per_level(test_grid,test_pft,:)
3078	!!$ WRITE(128,'(I8,20F20.10)') counter,gstot_component(test_grid,test_pft,:)
3079	!!$ WRITE(129,'(I8,20F20.10)') counter,gs_column_total(test_grid,test_pft)
3080	!!$ WRITE(130,'(I8,20F20.10)') counter,SUM(lai_per_level(test_grid,test_pft,:))
3081	!!$ IF(rveget(test_grid,test_pft) .LT. 1e19) &
3082	!!$ WRITE(131,'(I8,20ES20.10)') counter,rveget(test_grid,test_pft)
3083	!!$ IF(rveget_min(test_grid,test_pft) .LT. 1e19) &
3084	!!$ WRITE(132,'(I8,20ES20.10)') counter,rveget_min(test_grid,test_pft)
3085	!!$ WRITE(133,'(I8,20ES20.10)') counter,vbeta3(test_grid,test_pft)
3086	!!$ WRITE(134,'(I8,20ES20.10)') counter,vbeta3pot(test_grid,test_pft)
3087
3088	IF (printlev_loc>=3) WRITE (numout,*) ' diffuco_trans_co2 done '
3089
3090	END SUBROUTINE diffuco_trans_co2
3091
3092
3093	!! ================================================================================================================================
3094	!! SUBROUTINE : diffuco_comb
3095	!!
3096	!>\BRIEF This routine combines the previous partial beta
3097	!! coefficients and calculates the total alpha and complete beta coefficients.
3098	!!
3099	!! DESCRIPTION : Those integrated coefficients are used to calculate (in enerbil.f90) the total evapotranspiration
3100	!! from the grid-cell. \n
3101	!!
3102	!! In the case that air is more humid than surface, dew deposition can occur (negative latent heat flux).
3103	!! In this instance, for temperature above zero, all of the beta coefficients are set to 0, except for
3104	!! interception (vbeta2) and bare soil (vbeta4 with zero soil resistance). The amount of water that is
3105	!! intercepted by leaves is calculated based on the value of LAI of the surface. In the case of freezing
3106	!! temperatures, water is added to the snow reservoir, and so vbeta4 and vbeta2 are set to 0, and the
3107	!! total vbeta is set to 1.\n
3108	!!
3109	!! \latexonly
3110	!! \input{diffucocomb1.tex}
3111	!! \endlatexonly
3112	!!
3113	!! The beta and alpha coefficients are initially set to 1.
3114	!! \latexonly
3115	!! \input{diffucocomb2.tex}
3116	!! \endlatexonly
3117	!!
3118	!! If snow is lower than the critical value:
3119	!! \latexonly
3120	!! \input{diffucocomb3.tex}
3121	!! \endlatexonly
3122	!! If in the presence of dew:
3123	!! \latexonly
3124	!! \input{diffucocomb4.tex}
3125	!! \endlatexonly
3126	!!
3127	!! Determine where the water goes (soil, vegetation, or snow)
3128	!! when air moisture exceeds saturation.
3129	!! \latexonly
3130	!! \input{diffucocomb5.tex}
3131	!! \endlatexonly
3132	!!
3133	!! If it is not freezing dew is put into the interception reservoir and onto the bare soil. If it is freezing,
3134	!! water is put into the snow reservoir.
3135	!! Now modify vbetas where necessary: for soil and snow
3136	!! \latexonly
3137	!! \input{diffucocomb6.tex}
3138	!! \endlatexonly
3139	!!
3140	!! and for vegetation
3141	!! \latexonly
3142	!! \input{diffucocomb7.tex}
3143	!! \endlatexonly
3144	!!
3145	!! Then compute part of dew that can be intercepted by leafs.
3146	!!
3147	!! There will be no transpiration when air moisture is too high, under any circumstance
3148	!! \latexonly
3149	!! \input{diffucocomb8.tex}
3150	!! \endlatexonly
3151	!!
3152	!! There will also be no interception loss on bare soil, under any circumstance.
3153	!! \latexonly
3154	!! \input{diffucocomb9.tex}
3155	!! \endlatexonly
3156	!!
3157	!! The flowchart details the 'decision tree' which underlies the module.
3158	!!
3159	!! RECENT CHANGE(S): None
3160	!!
3161	!! MAIN OUTPUT VARIABLE(S): vbeta1, vbeta4, humrel, vbeta2, vbeta3, vbeta
3162	!!
3163	!! REFERENCE(S) :
3164	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
3165	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
3166	!! Circulation Model. Journal of Climate, 6, pp.248-273
3167	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
3168	!! sur le cycle de l'eau, PhD Thesis, available from:
3169	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
3170	!!
3171	!! FLOWCHART :
3172	!! \latexonly
3173	!! \includegraphics[scale=0.25]{diffuco_comb_flowchart.png}
3174	!! \endlatexonly
3175	!! \n
3176	!_ ================================================================================================================================
3177
3178	SUBROUTINE diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, &
3179	snow, veget, circ_class_biomass, circ_class_n,&
3180	tot_bare_soil, vbeta1, vbeta2, vbeta3 , vbeta4, vbeta, qsintmax, vbeta2sum, vbeta3sum)
3181
3182	! Ajout qsintmax dans les arguments de la routine Nathalie / le 13-03-2006
3183
3184	!! 0. Variable and parameter declaration
3185
3186	!! 0.1 Input variables
3187
3188	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (-)
3189	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
3190	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
3191	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Nortward Lowest level wind speed (m s^{-1})
3192	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
3193	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa)
3194	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
3195	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature (K)
3196	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lower air temperature (K)
3197	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
3198	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (fraction)
3199	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2})
3200	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tot_bare_soil!! Total evaporating bare soil fraction
3201	REAL(r_std),DIMENSION (kjpindex,nvm,ncirc,nparts,nelements), INTENT (in) :: circ_class_biomass !!
3202	REAL(r_std),DIMENSION (kjpindex,nvm,ncirc), INTENT (in) :: circ_class_n !!
3203	!! 0.2 Output variables
3204
3205	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient (-)
3206
3207	!! 0.3 Modified variables
3208
3209	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta1 !! Beta for sublimation (-)
3210	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta4 !! Beta for Bare soil evaporation (-)
3211	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Soil moisture stress (within range 0 to 1)
3212	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta2 !! Beta for interception loss (-)
3213	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3 !! Beta for Transpiration (-)
3214	REAL(r_std), DIMENSION(kjpindex), INTENT (inout) :: vbeta2sum, vbeta3sum
3215
3216	!! 0.4 Local variables
3217
3218	INTEGER(i_std) :: ji, jv
3219	REAL(r_std) :: zevtest, zsoil_moist, zrapp
3220	REAL(r_std), DIMENSION(kjpindex) :: qsatt
3221	LOGICAL, DIMENSION(kjpindex) :: toveg, tosnow
3222	REAL(r_std) :: coeff_dew_veg
3223	REAL(r_std), DIMENSION(kjpindex) :: vbeta2sum_temp,vbeta3sum_temp
3224	REAL(r_std),DIMENSION (kjpindex,nvm) :: lai !! Leaf area index (m^2 m^{-2})
3225	!_ ================================================================================================================================
3226
3227	!! \latexonly
3228	!! \input{diffucocomb1.tex}
3229	!! \endlatexonly
3230
3231	!! 0. Initialisation of vbeta2sum and vbeta3sum. They are needed in the
3232	!! calculations for the multi-layer energy budget.
3233	vbeta2sum(:) = zero
3234	vbeta3sum(:) = zero
3235
3236	!! 0.1 Calculate the LAI from the biomass.
3237	! This is done to prevent passing an incorrect
3238	! LAI around, since LAI is diagnostic.
3239	! We can't do this if we use an LAI map! So this breaks DO_STOMATE=N
3240	lai(:,ibare_sechiba) = zero
3241	DO jv = 1, nvm
3242	DO ji = 1, kjpindex
3243	lai(ji,jv) = cc_to_lai(circ_class_biomass(ji,jv,:,ileaf,icarbon),circ_class_n(ji,jv,:),jv)
3244	ENDDO
3245	ENDDO
3246
3247	!! 1. If dew is present:
3248
3249	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
3250
3251
3252	!! 1.1 Determine where the water goes
3253	!! Determine where the water goes (soil, vegetation, or snow)
3254	!! when air moisture exceeds saturation.
3255	!! \latexonly
3256	!! \input{diffucocomb5.tex}
3257	!! \endlatexonly
3258	toveg(:) = .FALSE.
3259	tosnow(:) = .FALSE.
3260	DO ji = 1, kjpindex
3261	IF ( qsatt(ji) .LT. qair(ji) ) THEN
3262	IF (temp_air(ji) .GT. tp_00) THEN
3263
3264	!! If it is not freezing dew is put into the
3265	!! interception reservoir and onto the bare soil.
3266	toveg(ji) = .TRUE.
3267	ELSE
3268
3269	!! If it is freezing water is put into the
3270	!! snow reservoir.
3271	tosnow(ji) = .TRUE.
3272	ENDIF
3273	ENDIF
3274	END DO
3275
3276
3277	!! 1.2 Now modify vbetas where necessary.
3278
3279	!! 1.2.1 Soil and snow
3280	!! \latexonly
3281	!! \input{diffucocomb6.tex}
3282	!! \endlatexonly
3283	DO ji = 1, kjpindex
3284	IF ( toveg(ji) ) THEN
3285	vbeta1(ji) = zero
3286	vbeta4(ji) = tot_bare_soil(ji)
3287	!! vbeta(ji) = vbeta4(ji)
3288	!! vbeta2sum(ji) = 0.0d0
3289	!! vbeta3sum(ji) = 0.0d0
3290	ENDIF
3291	IF ( tosnow(ji) ) THEN
3292	vbeta1(ji) = un
3293	vbeta4(ji) = zero
3294	!! vbeta(ji) = un
3295	!! vbeta2sum(ji) = 1.0d0
3296	!! vbeta3sum(ji) = 0.0d0
3297	ENDIF
3298	ENDDO
3299
3300	!! 1.2.2 Vegetation
3301	!! \latexonly
3302	!! \input{diffucocomb7.tex}
3303	!! \endlatexonly
3304	DO jv = 1, nvm
3305
3306	DO ji = 1, kjpindex
3307
3308	IF ( toveg(ji) ) THEN
3309	IF (qsintmax(ji,jv) .GT. min_sechiba) THEN
3310
3311	! Compute part of dew that can be intercepted by leafs.
3312	IF ( lai(ji,jv) .GT. min_sechiba) THEN
3313	IF (lai(ji,jv) .GT. 1.5) THEN
3314	coeff_dew_veg= &
3315	& dew_veg_poly_coeff(6)lai(ji,jv)*5 &
3316	& - dew_veg_poly_coeff(5)lai(ji,jv)*4 &
3317	& + dew_veg_poly_coeff(4)lai(ji,jv)*3 &
3318	& - dew_veg_poly_coeff(3)lai(ji,jv)*2 &
3319	& + dew_veg_poly_coeff(2)*lai(ji,jv) &
3320	& + dew_veg_poly_coeff(1)
3321	ELSE
3322	coeff_dew_veg=un
3323	ENDIF
3324	ELSE
3325	coeff_dew_veg=zero
3326	ENDIF
3327	IF (jv .EQ. 1) THEN
3328	! This line may not work with CWRR when frac_bare is distributed
3329	! among three soiltiles. Fortunately, qsintmax(ji,1)=0 (LAI=0 in
3330	! PFT1) so we never pass here
3331	vbeta2(ji,jv) = coeff_dew_veg*tot_bare_soil(ji)
3332	ELSE
3333	vbeta2(ji,jv) = coeff_dew_veg*veget(ji,jv)
3334	ENDIF
3335	ELSE
3336	vbeta2(ji,jv) = zero ! if qsintmax=0, vbeta2=0
3337	ENDIF
3338	!! vbeta(ji) = vbeta(ji) + vbeta2(ji,jv)
3339	!! vbeta4(ji) = 0.0d0
3340	!! vbeta2sum(ji) = vbeta(ji)
3341	!! vbeta3sum(ji) = 0.0d0
3342	ENDIF
3343	IF ( tosnow(ji) ) vbeta2(ji,jv) = zero
3344
3345	ENDDO
3346
3347	ENDDO
3348
3349	!! 1.2.3 Vegetation and transpiration
3350	!! There will be no transpiration when air moisture is too high, under any circumstance
3351	!! \latexonly
3352	!! \input{diffucocomb8.tex}
3353	!! \endlatexonly
3354	DO jv = 1, nvm
3355	DO ji = 1, kjpindex
3356	IF ( qsatt(ji) .LT. qair(ji) ) THEN
3357	vbeta3(ji,jv) = zero
3358	humrel(ji,jv) = zero
3359	ENDIF
3360	ENDDO
3361	ENDDO
3362
3363
3364	!! 1.2.4 Overrules 1.2.2
3365	!! There will also be no interception loss on bare soil, under any circumstance.
3366	!! \latexonly
3367	!! \input{diffucocomb9.tex}
3368	!! \endlatexonly
3369	DO ji = 1, kjpindex
3370	IF ( qsatt(ji) .LT. qair(ji) ) THEN
3371	vbeta2(ji,1) = zero
3372	ENDIF
3373	ENDDO
3374
3375	!! 2. Now calculate vbeta in all cases (the equality needs to hold for enerbil to be consistent)
3376	DO ji = 1, kjpindex
3377	vbeta(ji) = vbeta4(ji) + SUM(vbeta2(ji,:)) + SUM(vbeta3(ji,:))
3378
3379	vbeta2sum(ji) = SUM(vbeta2(ji,:))
3380	vbeta3sum(ji) = SUM(vbeta3(ji,:))
3381
3382
3383	IF (vbeta(ji) .LT. min_sechiba) THEN
3384	vbeta(ji) = zero
3385	vbeta4(ji) = zero
3386	vbeta2(ji,:)= zero
3387	vbeta3(ji,:)= zero
3388	END IF
3389	ENDDO
3390
3391	IF (printlev>=3) WRITE (numout,*) ' diffuco_comb done '
3392
3393	END SUBROUTINE diffuco_comb
3394
3395
3396	!! ================================================================================================================================
3397	!! SUBROUTINE : diffuco_raerod
3398	!!
3399	!>\BRIEF Computes the aerodynamic resistance, for cases in which the
3400	!! surface drag coefficient is provided by the coupled atmospheric model LMDZ and when the flag
3401	!! 'ldq_cdrag_from_gcm' is set to TRUE
3402	!!
3403	!! DESCRIPTION : Simply computes the aerodynamic resistance, for cases in which the
3404	!! surface drag coefficient is provided by the coupled atmospheric model LMDZ. If the surface drag coefficient
3405	!! is not provided by the LMDZ or signalled by the flag 'ldq_cdrag_from_gcm' set to FALSE, then the subroutine
3406	!! diffuco_aero is called instead of this one.
3407	!!
3408	!! Calculation of the aerodynamic resistance, for diganostic purposes. First calculate wind speed:
3409	!! \latexonly
3410	!! \input{diffucoaerod1.tex}
3411	!! \endlatexonly
3412	!!
3413	!! next calculate ::raero
3414	!! \latexonly
3415	!! \input{diffucoaerod2.tex}
3416	!! \endlatexonly
3417	!!
3418	!! RECENT CHANGE(S): None
3419	!!
3420	!! MAIN OUTPUT VARIABLE(S): ::raero
3421	!!
3422	!! REFERENCE(S) :
3423	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
3424	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
3425	!! Circulation Model. Journal of Climate, 6, pp.248-273
3426	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influence de l'irrigation
3427	!! sur le cycle de l'eau, PhD Thesis, available from:
3428	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
3429	!!
3430	!! FLOWCHART : None
3431	!! \n
3432	!_ ================================================================================================================================
3433
3434	SUBROUTINE diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
3435
3436	IMPLICIT NONE
3437
3438	!! 0. Variable and parameter declaration
3439
3440	!! 0.1 Input variables
3441
3442	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
3443	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind velocity (m s^{-1})
3444	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind velocity (m s^{-1})
3445	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
3446
3447	!! 0.2 Output variables
3448
3449	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: raero !! Aerodynamic resistance (s m^{-1})
3450
3451	!! 0.3 Modified variables
3452
3453	!! 0.4 Local variables
3454
3455	INTEGER(i_std) :: ji !! (-)
3456	REAL(r_std) :: speed !! (m s^{-1})
3457	!_ ================================================================================================================================
3458
3459	!! 1. Simple calculation of the aerodynamic resistance, for diganostic purposes.
3460
3461	DO ji=1,kjpindex
3462
3463	!! \latexonly
3464	!! \input{diffucoaerod1.tex}
3465	!! \endlatexonly
3466	speed = MAX(min_wind, wind(ji))
3467
3468	!! \latexonly
3469	!! \input{diffucoaerod2.tex}
3470	!! \endlatexonly
3471	raero(ji) = un / (q_cdrag(ji)*speed)
3472
3473	ENDDO
3474
3475	END SUBROUTINE diffuco_raerod
3476
3477
3478	FUNCTION Arrhenius (kjpindex,temp,ref_temp,energy_act) RESULT ( val_arrhenius )
3479	!! 0.1 Input variables
3480
3481	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
3482	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
3483	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
3484	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
3485
3486	!! 0.2 Result
3487
3488	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
3489	!! a Arrhenius function (-)
3490
3491	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:))))
3492	END FUNCTION Arrhenius
3493
3494	FUNCTION Arrhenius_modified_1d (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius )
3495	!! 0.1 Input variables
3496
3497	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
3498	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
3499	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
3500	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
3501	REAL(r_std),INTENT(in) :: energy_deact !! Deactivation Energy (J mol-1)
3502	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: entropy !! Entropy term (J K-1 mol-1)
3503
3504	!! 0.2 Result
3505
3506	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
3507	!! a Arrhenius function (-)
3508
3509	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:)))) &
3510	* (1. + EXP( (ref_temp * entropy(:) - energy_deact) / (ref_temp * RR ))) &
3511	/ (1. + EXP( (temp(:) * entropy(:) - energy_deact) / ( RR*temp(:))))
3512
3513	END FUNCTION Arrhenius_modified_1d
3514
3515	FUNCTION Arrhenius_modified_0d (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius )
3516	!! 0.1 Input variables
3517
3518	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
3519	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
3520	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
3521	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
3522	REAL(r_std),INTENT(in) :: energy_deact !! Deactivation Energy (J mol-1)
3523	REAL(r_std),INTENT(in) :: entropy !! Entropy term (J K-1 mol-1)
3524
3525	!! 0.2 Result
3526
3527	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
3528	!! a Arrhenius function (-)
3529
3530	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:)))) &
3531	* (1. + EXP( (ref_temp * entropy - energy_deact) / (ref_temp * RR ))) &
3532	/ (1. + EXP( (temp(:) * entropy - energy_deact) / ( RR*temp(:))))
3533
3534	END FUNCTION Arrhenius_modified_0d
3535
3536
3537	!! ================================================================================================================================
3538	!! SUBROUTINE : isotope_c13
3539	!!
3540	!>\BRIEF Calculate the fractionation of C13 during photosynthesis\n
3541	!!
3542	!! DESCRIPTION :
3543	!!
3544	!! RECENT CHANGE(S): None
3545	!!
3546	!! MAIN OUTPUT VARIABLE(S): :: delta_c13_assim
3547	!!
3548	!! REFERENCE(S) :
3549	!! - Farquhar, G. D., Ehleringer, J. R., & Hubick, K. T. (1989). CARBON ISOTOPE DISCRIMINATION AND PHOTOSYNTHESIS.
3550	!! Annual Review of Plant Physiology and Plant Molecular Biology, 40, 503-537. doi: 10.1146/annurev.arplant.40.1.50
3551	!!
3552	!! FLOWCHART : None
3553	!_ ================================================================================================================================
3554
3555	SUBROUTINE isotope_c13(kjpindex, index_assi, leaf_ci, ccanopy, lai_per_level, &
3556	assimtot, delta_c13_assim, leaf_ci_out, nvm, &
3557	nlevels_tot, assimi_lev)
3558
3559	IMPLICIT NONE
3560
3561	!! 0. Variables and parameter declaration
3562
3563	!! 0.1 Input variables
3564	INTEGER(i_std), INTENT (in) :: kjpindex !! number of land pixels
3565	INTEGER(i_std), INTENT (in) :: nvm !! number of pfts
3566	INTEGER(i_std), INTENT (in) :: nlevels_tot !! number of lai layers in canopy
3567	INTEGER(i_std), DIMENSION(:), INTENT (in) :: index_assi !! indices (unitless)
3568	REAL(r_std), DIMENSION (:,:,:), INTENT (in) :: leaf_ci !! intercellular CO2 concentration (ppm)
3569	REAL(r_std),DIMENSION (:), INTENT (in) :: ccanopy !! CO2 concentration inside the canopy (ppm)
3570	REAL(r_std), DIMENSION(:,:,:), INTENT (in) :: lai_per_level !! This is the LAI per vertical level
3571	!! @tex $(m^{2} m^{-2})$ @endtex
3572	!+++CHECK+++
3573	! Looks like assimi is a local variable in trans_co2 and that
3574	! its values are overwritten for every new layer. If this is
3575	! indeed the case, it should not be used outside trans_co2
3576	! REAL(r_std), DIMENSION(:,:), INTENT (in) :: assimi !! assimilation (at a specific LAI level)
3577	! !! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
3578	!+++++++++++
3579	REAL(r_std), DIMENSION(:,:), INTENT (in) :: assimtot !! total assimilation
3580	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
3581	REAL(r_std), DIMENSION(:,:,:), INTENT (in) :: assimi_lev !! assimilation (at all LAI levels)
3582	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
3583
3584	!! 0.2 Output variables
3585	REAL(r_std),DIMENSION (:,:), INTENT (out) :: delta_c13_assim !! C13 concentration in delta notation
3586	!! @tex $ permille $ @endtex (per thousand)
3587	REAL(r_std),DIMENSION (:,:), INTENT (out) :: leaf_ci_out !! Ci @tex $ ??unit?? $ @endtex
3588
3589	!! 0.3 Modified variables
3590
3591	!! 0.4 Local variables
3592	INTEGER(i_std) :: inia,iainia !! counter/indices (unitless)
3593	INTEGER(i_std) :: ji,jv,jl !! Indices (-)
3594
3595
3596	!_ ================================================================================================================================
3597
3598	!! Initialize
3599	delta_c13_assim(:,:) = zero
3600	leaf_ci_out(:,:) = zero
3601	!printlev_loc = get_printlev('c13')
3602
3603	!! Calculate delta_C13
3604	DO ji=1,kjpindex
3605
3606	! Skip the whole pixel of none of the PFTs did photosynthesis
3607	IF (SUM(SUM(assimi_lev(ji,:,:),2),1) .LE. zero) THEN
3608	! There was no photsynthesis at this pixel, no reason to do the
3609	! calculation
3610	CYCLE
3611	ENDIF
3612
3613	DO jv = 1,nvm
3614
3615	! Skip the PFT is it didn't do photosynthesis
3616	IF (SUM(assimi_lev(ji,jv,:)) .LE. zero) THEN
3617	! There was no photsynthesis at this pixel, no reason to do the
3618	! calculation
3619	CYCLE
3620	ENDIF
3621
3622	! There was photosynthesis at this PFT, calculate the C13 fractionation
3623	DO jl = 1,nlevels_tot
3624
3625	IF (lai_per_level(ji,jv,jl) .GT. zero) THEN
3626
3627	IF (printlev_loc>=4) THEN
3628	WRITE(numout,*) 'leaf_ci, ',leaf_ci(ji,jv,jl)
3629	WRITE(numout,*) 'assimi_lev, ',assimi_lev(ji,jv,jl)
3630	WRITE(numout,*) 'assimtot - isotope - before calc, ',assimtot(ji,jv)
3631	WRITE(numout,*) 'lai_per_level, ',lai_per_level(ji,jv,jl)
3632	WRITE(numout,*) 'day, ',day
3633	WRITE(numout,*) 'pft, ',jv
3634	ENDIF
3635
3636
3637	! Simple version of Farquhar model
3638	! delta_c13_assim(:,jv) = (-8 - c13_a - (c13_b - c13_a) * cimean(:,jv) / ccanopy(:))
3639	delta_c13_assim(ji,jv) = delta_c13_assim(ji,jv) + (- c13_a - (c13_b - c13_a)* &
3640	leaf_ci(ji,jv,jl)/ccanopy(ji)) * &
3641	assimi_lev(ji,jv,jl) * lai_per_level(ji,jv,jl)
3642	leaf_ci_out(ji,jv) = leaf_ci_out(ji,jv) + &
3643	leaf_ci(ji,jv,jl) * &
3644	assimi_lev(ji,jv,jl) * lai_per_level(ji,jv,jl)
3645
3646	ENDIF
3647
3648	! Simple version of Farquhar model with photorespiratioin, gi and gamma*
3649	! This code copied from the work by Thomas Eglin and Thomas Launois where it
3650	! was already commented out. When moving it to ORCHIDEE-CAN is has not been
3651	! tested with ORCHIDEE-CAN. It is not clear why it was commented out.
3652	! was commented out
3653	!!$ IF (assimi(ji,jv) .GT. 0) THEN
3654	!!$ delta_c13_assim(ji,jv) = delta_c13_assim(ji,jv) + (&
3655	!!$ - c13_ab * (ccanopy(ji)-leaf_cs(ji,jv,jl))/(ccanopy(ji)) &
3656	!!$ - c13_a * (leaf_cs(ji,jv,jl)-leaf_ci(ji,jv,jl))/(ccanopy(ji)) &
3657	!!$ - (C13_es + c13_al) * (leaf_ci(ji,jv,jl)-leaf_cc(ji,jv,jl))/(ccanopy(ji)) &
3658	!!$ - (c13_b+2.0) * leaf_cc(ji,jv,jl)/(ccanopy(ji)) &
3659	!!$ - (C13_e * (Rd(ji)/(Rd(ji) + assimi(ji,jv))) + C13_f * CP(ji))/(ccanopy(ji))&
3660	!!$ ) * assimi(ji,jv) * (laitab(jl+1)-laitab(jl))
3661	!!$ ENDIF
3662
3663
3664	!!$ (Rd(ji) + assimi(ji,jv)), Rd(ji),assimi(ji,jv)
3665
3666	! This code copied from the work by Thomas Eglin and Thomas Launois where it
3667	! was already commented out. When moving it to ORCHIDEE-CAN is has not been
3668	! tested with ORCHIDEE-CAN. It is not clear why it was commented out.
3669	! was commented out
3670	! Modification M; Bordigoni 07/2008 : New calculation of c3-plants fractionation
3671	! according to Lloyd and Farquhar_Oecologia_1994
3672	!!$ delta_c13_assim(ji,jv) = delta_c13_assim(ji,jv) +(- (&
3673	!!$ 4.4*(1-(leaf_ci(ji,jv,jl) / ccanopy(ji) )+0.025)+ &
3674	!!$ 0.075*(1.1+0.7)+ &
3675	!!$ 27.5*(( leaf_ci(ji,jv,jl) / ccanopy(ji))-0.1)- &
3676	!!$ (110.1+8(1.541.05((temp_air(ji)-273.16)+2.5)))/ ccanopy(ji)) )* &
3677	!!$ assimi(ji,jv) * (laitab(jl+1)-laitab(jl))
3678
3679	ENDDO ! photosynthesis layers
3680
3681	IF (assimtot(ji,jv) .GT. threshold_c13_assim) THEN
3682
3683	delta_c13_assim(ji,jv) = delta_c13_assim(ji,jv)/assimtot(ji,jv)
3684	leaf_ci_out(ji,jv) = leaf_ci_out(ji,jv)/assimtot(ji,jv)
3685
3686	ELSE
3687
3688	delta_c13_assim(ji,jv) = zero
3689	leaf_ci_out(ji,jv) = zero
3690
3691	ENDIF
3692
3693	ENDDO ! # PFTs
3694
3695	ENDDO ! #kjpindex
3696
3697	END SUBROUTINE isotope_c13
3698
3699	END MODULE diffuco

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

Download in other formats: