1 | ! ================================================================================================================================ |
---|
2 | ! MODULE : diffuco |
---|
3 | ! |
---|
4 | ! CONTACT : orchidee-help _at_ listes.ipsl.fr |
---|
5 | ! |
---|
6 | ! LICENCE : IPSL (2006) |
---|
7 | ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
8 | ! |
---|
9 | !>\BRIEF 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 constantes |
---|
54 | USE constantes_soil |
---|
55 | USE qsat_moisture |
---|
56 | USE sechiba_io_p |
---|
57 | USE ioipsl |
---|
58 | USE pft_parameters |
---|
59 | USE grid |
---|
60 | USE time, ONLY : one_day, dt_sechiba |
---|
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 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: wind !! Wind module (m s^{-1}) |
---|
78 | !$OMP THREADPRIVATE(wind) |
---|
79 | |
---|
80 | CONTAINS |
---|
81 | |
---|
82 | |
---|
83 | !! ============================================================================================================================= |
---|
84 | !! SUBROUTINE: diffuco_initialize |
---|
85 | !! |
---|
86 | !>\BRIEF Allocate module variables, read from restart file or initialize with default values |
---|
87 | !! |
---|
88 | !! DESCRIPTION: Allocate module variables, read from restart file or initialize with default values. |
---|
89 | !! Call chemistry_initialize for initialization of variables needed for the calculations of BVOCs. |
---|
90 | !! |
---|
91 | !! RECENT CHANGE(S): None |
---|
92 | !! |
---|
93 | !! REFERENCE(S): None |
---|
94 | !! |
---|
95 | !! FLOWCHART: None |
---|
96 | !! \n |
---|
97 | !_ ============================================================================================================================== |
---|
98 | SUBROUTINE diffuco_initialize (kjit, kjpindex, index, & |
---|
99 | rest_id, lalo, neighbours, resolution, & |
---|
100 | rstruct, q_cdrag) |
---|
101 | |
---|
102 | !! 0. Variable and parameter declaration |
---|
103 | !! 0.1 Input variables |
---|
104 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-) |
---|
105 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
106 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-) |
---|
107 | INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (-) |
---|
108 | REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates |
---|
109 | INTEGER(i_std),DIMENSION (kjpindex,NbNeighb),INTENT (in):: neighbours !! Vector of neighbours for each |
---|
110 | REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y |
---|
111 | |
---|
112 | !! 0.2 Output variables |
---|
113 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! Structural resistance for the vegetation |
---|
114 | |
---|
115 | !! 0.3 Modified variables |
---|
116 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
117 | |
---|
118 | !! 0.4 Local variables |
---|
119 | INTEGER :: ilai |
---|
120 | INTEGER :: jv |
---|
121 | INTEGER :: ier |
---|
122 | CHARACTER(LEN=4) :: laistring |
---|
123 | CHARACTER(LEN=80) :: var_name |
---|
124 | !_ ================================================================================================================================ |
---|
125 | |
---|
126 | !! 1. Define flag ldq_cdrag_from_gcm. This flag determines if the cdrag should be taken from the GCM or be calculated. |
---|
127 | !! The default value is true if the q_cdrag variables was already initialized. This is the case when coupled to the LMDZ. |
---|
128 | |
---|
129 | !Config Key = CDRAG_FROM_GCM |
---|
130 | !Config Desc = Keep cdrag coefficient from gcm. |
---|
131 | !Config If = OK_SECHIBA |
---|
132 | !Config Def = y |
---|
133 | !Config Help = Set to .TRUE. if you want q_cdrag coming from GCM (if q_cdrag on initialization is non zero). |
---|
134 | !Config Keep cdrag coefficient from gcm for latent and sensible heat fluxes. |
---|
135 | !Config Units = [FLAG] |
---|
136 | IF ( ABS(MAXVAL(q_cdrag)) .LE. EPSILON(q_cdrag)) THEN |
---|
137 | ldq_cdrag_from_gcm = .FALSE. |
---|
138 | ELSE |
---|
139 | ldq_cdrag_from_gcm = .TRUE. |
---|
140 | ENDIF |
---|
141 | CALL getin_p('CDRAG_from_GCM', ldq_cdrag_from_gcm) |
---|
142 | IF (printlev>=2) WRITE(numout,*) "ldq_cdrag_from_gcm = ",ldq_cdrag_from_gcm |
---|
143 | |
---|
144 | !! 2. Allocate module variables |
---|
145 | ALLOCATE (wind(kjpindex),stat=ier) |
---|
146 | IF (ier /= 0) CALL ipslerr_p(3,'diffuco_initialize','Problem in allocate of variable wind','','') |
---|
147 | |
---|
148 | !! 3. Read variables from restart file |
---|
149 | IF (printlev>=3) WRITE (numout,*) 'Read DIFFUCO variables from restart file' |
---|
150 | |
---|
151 | CALL ioconf_setatt_p('UNITS', 's/m') |
---|
152 | CALL ioconf_setatt_p('LONG_NAME','Structural resistance') |
---|
153 | CALL restget_p (rest_id, 'rstruct', nbp_glo, nvm, 1, kjit, .TRUE., rstruct, "gather", nbp_glo, index_g) |
---|
154 | IF ( ALL(rstruct(:,:) == val_exp) ) THEN |
---|
155 | DO jv = 1, nvm |
---|
156 | rstruct(:,jv) = rstruct_const(jv) |
---|
157 | ENDDO |
---|
158 | ENDIF |
---|
159 | |
---|
160 | !! 4. Initialize chemistry module |
---|
161 | IF (printlev>=3) WRITE(numout,*) "ok_bvoc:",ok_bvoc |
---|
162 | IF ( ok_bvoc ) CALL chemistry_initialize(kjpindex, lalo, neighbours, resolution) |
---|
163 | |
---|
164 | END SUBROUTINE diffuco_initialize |
---|
165 | |
---|
166 | |
---|
167 | |
---|
168 | !! ================================================================================================================================ |
---|
169 | !! SUBROUTINE : diffuco_main |
---|
170 | !! |
---|
171 | !>\BRIEF The root subroutine for the module, which calls all other required |
---|
172 | !! subroutines. |
---|
173 | !! |
---|
174 | !! DESCRIPTION : |
---|
175 | |
---|
176 | !! This is the main subroutine for the module. |
---|
177 | !! First it calculates the surface drag coefficient (via a call to diffuco_aero), using available parameters to determine |
---|
178 | !! the stability of air in the surface layer by calculating the Richardson Nubmber. If a value for the |
---|
179 | !! surface drag coefficient is passed down from the atmospheric model and and if the flag 'ldq_cdrag_from_gcm' |
---|
180 | !! is set to TRUE, then the subroutine diffuco_aerod is called instead. This calculates the aerodynamic coefficient. \n |
---|
181 | !! |
---|
182 | !! Following this, an estimation of the saturated humidity at the surface is made (via a call |
---|
183 | !! to qsatcalc in the module qsat_moisture). Following this the beta coefficients for sublimation (via |
---|
184 | !! diffuco_snow), interception (diffuco_inter), bare soil (diffuco_bare), and transpiration (via |
---|
185 | !! diffuco_trans_co2) are calculated in sequence. Finally |
---|
186 | !! the alpha and beta coefficients are combined (diffuco_comb). \n |
---|
187 | !! |
---|
188 | !! The surface drag coefficient is calculated for use within the module enerbil. It is required to to |
---|
189 | !! calculate the aerodynamic coefficient for every flux. \n |
---|
190 | !! |
---|
191 | !! The various beta coefficients are used within the module enerbil for modifying the rate of evaporation, |
---|
192 | !! as appropriate for the surface. As explained in Chapter 2 of Guimberteau (2010), that module (enerbil) |
---|
193 | !! calculates the rate of evaporation essentially according to the expression $E = /beta E_{pot}$, where |
---|
194 | !! E is the total evaporation and $E_{pot}$ the evaporation potential. If $\beta = 1$, there would be |
---|
195 | !! essentially no resistance to evaporation, whereas should $\beta = 0$, there would be no evaporation and |
---|
196 | !! the surface layer would be subject to some very stong hydrological stress. \n |
---|
197 | !! |
---|
198 | !! The following processes are calculated: |
---|
199 | !! - call diffuco_aero for aerodynamic transfer coeficient |
---|
200 | !! - call diffuco_snow for partial beta coefficient: sublimation |
---|
201 | !! - call diffuco_inter for partial beta coefficient: interception for each type of vegetation |
---|
202 | !! - call diffuco_bare for partial beta coefficient: bare soil |
---|
203 | !! - call diffuco_trans_co2 for partial beta coefficient: transpiration for each type of vegetation, using Farqhar's formula |
---|
204 | !! - call diffuco_comb for alpha and beta coefficient |
---|
205 | !! - call chemistry_bvoc for alpha and beta coefficients for biogenic emissions |
---|
206 | !! |
---|
207 | !! RECENT CHANGE(S): None |
---|
208 | !! |
---|
209 | !! MAIN OUTPUT VARIABLE(S): humrel, q_cdrag, vbeta, vbeta1, vbeta4, |
---|
210 | !! vbeta2, vbeta3, rveget, cimean |
---|
211 | !! |
---|
212 | !! REFERENCE(S) : |
---|
213 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
214 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
215 | !! Circulation Model. Journal of Climate, 6, pp.248-273. |
---|
216 | !! - de Rosnay, P, 1999. Représentation des interactions sol-plante-atmosphÚre dans le modÚle de circulation générale |
---|
217 | !! du LMD, 1999. PhD Thesis, Université Paris 6, available (25/01/12): |
---|
218 | !! http://www.ecmwf.int/staff/patricia_de_rosnay/publications.html#8 |
---|
219 | !! - Ducharne, A, 1997. Le cycle de l'eau: modélisation de l'hydrologie continentale, étude de ses interactions avec |
---|
220 | !! le climat, PhD Thesis, Université Paris 6 |
---|
221 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation |
---|
222 | !! sur le cycle de l'eau, PhD Thesis, available (25/01/12): |
---|
223 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf |
---|
224 | !! - LathiÚre, J, 2005. Evolution des émissions de composés organiques et azotés par la biosphÚre continentale dans le |
---|
225 | !! modÚle LMDz-INCA-ORCHIDEE, Université Paris 6 |
---|
226 | !! |
---|
227 | !! FLOWCHART : |
---|
228 | !! \latexonly |
---|
229 | !! \includegraphics[scale=0.5]{diffuco_main_flowchart.png} |
---|
230 | !! \endlatexonly |
---|
231 | !! \n |
---|
232 | !_ ================================================================================================================================ |
---|
233 | |
---|
234 | SUBROUTINE diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, & |
---|
235 | & zlev, z0m, z0h, roughheight, temp_sol, temp_air, temp_growth, rau, q_cdrag, qsurf, qair, pb, & |
---|
236 | & evap_bare_lim, evap_bare_lim_ns, evapot, evapot_corr, snow, flood_frac, flood_res, & |
---|
237 | & frac_nobio, snow_nobio, totfrac_nobio, & |
---|
238 | & swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, & |
---|
239 | & vbeta , vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, gsmean, rveget, rstruct, cimean, gpp, & |
---|
240 | & lalo, neighbours, resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, frac_snow_veg, frac_snow_nobio, & |
---|
241 | & hist_id, hist2_id) |
---|
242 | |
---|
243 | !! 0. Variable and parameter declaration |
---|
244 | |
---|
245 | !! 0.1 Input variables |
---|
246 | |
---|
247 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-) |
---|
248 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
249 | INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (-) |
---|
250 | INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier (-) |
---|
251 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-) |
---|
252 | INTEGER(i_std),DIMENSION (kjpindex*(nlai+1)), INTENT (in) :: indexlai !! Indeces of the points on the 3D map |
---|
253 | INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map (-) |
---|
254 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1}) |
---|
255 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1}) |
---|
256 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m) |
---|
257 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0m !! Surface roughness Length for momentum (m) |
---|
258 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0h !! Surface roughness Length for heat (m) |
---|
259 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective height for roughness (m) |
---|
260 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature (K) |
---|
261 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature (K) |
---|
262 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_growth !! Growth temperature (°C) - Is equal to t2m_month |
---|
263 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3}) |
---|
264 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Near surface air specific humidity (kg kg^{-1}) |
---|
265 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level air specific humidity (kg kg^{-1}) |
---|
266 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg) |
---|
267 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_frac !! Fraction of floodplains |
---|
268 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_res !! Reservoir in floodplains (estimation to avoid over-evaporation) |
---|
269 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface level pressure (hPa) |
---|
270 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evap_bare_lim !! Limit to the bare soil evaporation when the |
---|
271 | !! 11-layer hydrology is used (-) |
---|
272 | REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout) :: evap_bare_lim_ns !! Limit to the bare soil evaporation when the |
---|
273 | !! 11-layer hydrology is used (-) |
---|
274 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm day^{-1}) |
---|
275 | !! NdN This variable does not seem to be used at |
---|
276 | !! all in diffuco |
---|
277 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation |
---|
278 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice,lakes,cities,... (-) |
---|
279 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice,lakes,cities,... (kg m^{-2}) |
---|
280 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+cities+... (-) |
---|
281 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2}) |
---|
282 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux (W m^{-2}) |
---|
283 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: coszang !! Cosine of the solar zenith angle (unitless) |
---|
284 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration inside the canopy (ppm) |
---|
285 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (-) |
---|
286 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty) |
---|
287 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2}) |
---|
288 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2}) |
---|
289 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception |
---|
290 | !! (kg m^{-2}) |
---|
291 | REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param !! referenced vcmax, nue, leaf N for photosynthesis |
---|
292 | REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates |
---|
293 | INTEGER(i_std),DIMENSION (kjpindex,NbNeighb),INTENT (in):: neighbours !! Vector of neighbours for each |
---|
294 | !! grid point (1=N, 2=E, 3=S, 4=W) |
---|
295 | REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y |
---|
296 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ptnlev1 !! 1st level of soil temperature (Kelvin) |
---|
297 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation expressed in mm/tstep |
---|
298 | REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT (in) :: frac_age !! Age efficiency for isoprene emissions (from STOMATE) |
---|
299 | REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction |
---|
300 | REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area |
---|
301 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in):: frac_snow_nobio !! Snow cover fraction on non-vegeted area |
---|
302 | |
---|
303 | !! 0.2 Output variables |
---|
304 | |
---|
305 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient (-) |
---|
306 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation (-) |
---|
307 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation (-) |
---|
308 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta5 !! Beta for floodplains |
---|
309 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gsmean !! Mean stomatal conductance to CO2 (mol m-2 s-1) |
---|
310 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss (-) |
---|
311 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for transpiration (-) |
---|
312 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3pot !! Beta for potential transpiration |
---|
313 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Stomatal resistance for the whole canopy (s m^{-1}) |
---|
314 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! Structural resistance for the vegetation |
---|
315 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! Mean leaf Ci (ppm) |
---|
316 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out) :: gpp !! Assimilation ((gC m^{-2} dt_sechiba^{-1}), total area) |
---|
317 | |
---|
318 | !! 0.3 Modified variables |
---|
319 | |
---|
320 | REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (inout) :: humrel !! Soil moisture stress (within range 0 to 1) |
---|
321 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Surface drag coefficient (-) |
---|
322 | |
---|
323 | !! 0.4 Local variables |
---|
324 | |
---|
325 | REAL(r_std),DIMENSION (kjpindex,nvm) :: vbeta23 !! Beta for fraction of wetted foliage that will |
---|
326 | !! transpire once intercepted water has evaporated (-) |
---|
327 | REAL(r_std),DIMENSION (kjpindex) :: raero !! Aerodynamic resistance (s m^{-1}) |
---|
328 | INTEGER(i_std) :: ilai |
---|
329 | CHARACTER(LEN=4) :: laistring |
---|
330 | CHARACTER(LEN=80) :: var_name !! To store variables names for I/O |
---|
331 | REAL(r_std),DIMENSION(kjpindex) :: qsatt !! Surface saturated humidity (kg kg^{-1}) |
---|
332 | REAL(r_std),DIMENSION (kjpindex,nvm) :: cim !! Intercellular CO2 over nlai |
---|
333 | |
---|
334 | !_ ================================================================================================================================ |
---|
335 | |
---|
336 | wind(:) = SQRT (u(:)*u(:) + v(:)*v(:)) |
---|
337 | |
---|
338 | !! 1. Calculate the different coefficients |
---|
339 | |
---|
340 | IF (.NOT.ldq_cdrag_from_gcm) THEN |
---|
341 | ! Case 1a) |
---|
342 | CALL diffuco_aero (kjpindex, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, & |
---|
343 | qsurf, qair, snow, q_cdrag) |
---|
344 | ENDIF |
---|
345 | |
---|
346 | ! Case 1b) |
---|
347 | CALL diffuco_raerod (kjpindex, u, v, q_cdrag, raero) |
---|
348 | |
---|
349 | !! 2. Make an estimation of the saturated humidity at the surface |
---|
350 | |
---|
351 | CALL qsatcalc (kjpindex, temp_sol, pb, qsatt) |
---|
352 | |
---|
353 | !! 3. Calculate the beta coefficient for sublimation |
---|
354 | |
---|
355 | CALL diffuco_snow (kjpindex, qair, qsatt, rau, u, v, q_cdrag, & |
---|
356 | snow, frac_nobio, totfrac_nobio, snow_nobio, frac_snow_veg, frac_snow_nobio, & |
---|
357 | vbeta1) |
---|
358 | |
---|
359 | |
---|
360 | CALL diffuco_flood (kjpindex, qair, qsatt, rau, u, v, q_cdrag, evapot, evapot_corr, & |
---|
361 | & flood_frac, flood_res, vbeta5) |
---|
362 | |
---|
363 | !! 4. Calculate the beta coefficient for interception |
---|
364 | |
---|
365 | CALL diffuco_inter (kjpindex, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, & |
---|
366 | & qsintveg, qsintmax, rstruct, vbeta2, vbeta23) |
---|
367 | |
---|
368 | |
---|
369 | !! 5. Calculate the beta coefficient for transpiration |
---|
370 | |
---|
371 | CALL diffuco_trans_co2 (kjpindex, swdown, pb, qsurf, qair, temp_air, temp_growth, rau, u, v, q_cdrag, humrel, & |
---|
372 | assim_param, ccanopy, & |
---|
373 | veget, veget_max, lai, qsintveg, qsintmax, vbeta3, vbeta3pot, & |
---|
374 | rveget, rstruct, cimean, gsmean, gpp, & |
---|
375 | vbeta23, hist_id, indexveg, indexlai, index, kjit, cim) |
---|
376 | |
---|
377 | ! |
---|
378 | !biogenic emissions |
---|
379 | ! |
---|
380 | IF ( ok_bvoc ) THEN |
---|
381 | CALL chemistry_bvoc (kjpindex, swdown, coszang, temp_air, & |
---|
382 | temp_sol, ptnlev1, precip_rain, humrel, veget_max, & |
---|
383 | lai, frac_age, lalo, ccanopy, cim, wind, snow, & |
---|
384 | veget, hist_id, hist2_id, kjit, index, & |
---|
385 | indexlai, indexveg) |
---|
386 | ENDIF |
---|
387 | ! |
---|
388 | ! combination of coefficient : alpha and beta coefficient |
---|
389 | ! beta coefficient for bare soil |
---|
390 | ! |
---|
391 | |
---|
392 | CALL diffuco_bare (kjpindex, tot_bare_soil, veget_max, evap_bare_lim, evap_bare_lim_ns, vbeta2, vbeta3, vbeta4) |
---|
393 | |
---|
394 | !! 6. Combine the alpha and beta coefficients |
---|
395 | |
---|
396 | ! Ajout qsintmax dans les arguments de la routine.... Nathalie / le 13-03-2006 |
---|
397 | CALL diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, snow, & |
---|
398 | & veget, lai, tot_bare_soil, vbeta1, vbeta2, vbeta3, vbeta4, & |
---|
399 | & evap_bare_lim, evap_bare_lim_ns, veget_max, vbeta, qsintmax) |
---|
400 | |
---|
401 | CALL xios_orchidee_send_field("q_cdrag",q_cdrag) |
---|
402 | CALL xios_orchidee_send_field("raero",raero) |
---|
403 | CALL xios_orchidee_send_field("wind",wind) |
---|
404 | CALL xios_orchidee_send_field("qsatt",qsatt) |
---|
405 | CALL xios_orchidee_send_field("coszang",coszang) |
---|
406 | CALL xios_orchidee_send_field('cim', cim) |
---|
407 | |
---|
408 | IF ( .NOT. almaoutput ) THEN |
---|
409 | CALL histwrite_p(hist_id, 'raero', kjit, raero, kjpindex, index) |
---|
410 | CALL histwrite_p(hist_id, 'cdrag', kjit, q_cdrag, kjpindex, index) |
---|
411 | CALL histwrite_p(hist_id, 'Wind', kjit, wind, kjpindex, index) |
---|
412 | CALL histwrite_p(hist_id, 'qsatt', kjit, qsatt, kjpindex, index) |
---|
413 | CALL histwrite_p(hist_id, 'cim', kjit, cim, kjpindex*nvm, indexveg) |
---|
414 | |
---|
415 | IF ( hist2_id > 0 ) THEN |
---|
416 | CALL histwrite_p(hist2_id, 'raero', kjit, raero, kjpindex, index) |
---|
417 | CALL histwrite_p(hist2_id, 'cdrag', kjit, q_cdrag, kjpindex, index) |
---|
418 | CALL histwrite_p(hist2_id, 'Wind', kjit, wind, kjpindex, index) |
---|
419 | CALL histwrite_p(hist2_id, 'qsatt', kjit, qsatt, kjpindex, index) |
---|
420 | ENDIF |
---|
421 | ELSE |
---|
422 | CALL histwrite_p(hist_id, 'cim', kjit, cim, kjpindex*nvm, indexveg) |
---|
423 | ENDIF |
---|
424 | |
---|
425 | IF (printlev>=3) WRITE (numout,*) ' diffuco_main done ' |
---|
426 | |
---|
427 | END SUBROUTINE diffuco_main |
---|
428 | |
---|
429 | !! ============================================================================================================================= |
---|
430 | !! SUBROUTINE: diffuco_finalize |
---|
431 | !! |
---|
432 | !>\BRIEF Write to restart file |
---|
433 | !! |
---|
434 | !! DESCRIPTION: This subroutine writes the module variables and variables calculated in diffuco |
---|
435 | !! to restart file |
---|
436 | !! |
---|
437 | !! RECENT CHANGE(S): None |
---|
438 | !! REFERENCE(S): None |
---|
439 | !! FLOWCHART: None |
---|
440 | !! \n |
---|
441 | !_ ============================================================================================================================== |
---|
442 | SUBROUTINE diffuco_finalize (kjit, kjpindex, rest_id, rstruct ) |
---|
443 | |
---|
444 | !! 0. Variable and parameter declaration |
---|
445 | !! 0.1 Input variables |
---|
446 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-) |
---|
447 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
448 | INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (-) |
---|
449 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! Structural resistance for the vegetation |
---|
450 | |
---|
451 | !! 0.4 Local variables |
---|
452 | INTEGER :: ilai |
---|
453 | CHARACTER(LEN=4) :: laistring |
---|
454 | CHARACTER(LEN=80) :: var_name |
---|
455 | |
---|
456 | !_ ================================================================================================================================ |
---|
457 | |
---|
458 | !! 1. Prepare the restart file for the next simulation |
---|
459 | IF (printlev>=3) WRITE (numout,*) 'Complete restart file with DIFFUCO variables ' |
---|
460 | |
---|
461 | CALL restput_p (rest_id, 'rstruct', nbp_glo, nvm, 1, kjit, rstruct, 'scatter', nbp_glo, index_g) |
---|
462 | |
---|
463 | END SUBROUTINE diffuco_finalize |
---|
464 | |
---|
465 | |
---|
466 | !! ================================================================================================================================ |
---|
467 | !! SUBROUTINE : diffuco_clear |
---|
468 | !! |
---|
469 | !>\BRIEF Housekeeping module to deallocate the variables |
---|
470 | !! rstruct and raero |
---|
471 | !! |
---|
472 | !! DESCRIPTION : Housekeeping module to deallocate the variables |
---|
473 | !! rstruct and raero |
---|
474 | !! |
---|
475 | !! RECENT CHANGE(S) : None |
---|
476 | !! |
---|
477 | !! MAIN OUTPUT VARIABLE(S) : None |
---|
478 | !! |
---|
479 | !! REFERENCE(S) : None |
---|
480 | !! |
---|
481 | !! FLOWCHART : None |
---|
482 | !! \n |
---|
483 | !_ ================================================================================================================================ |
---|
484 | |
---|
485 | SUBROUTINE diffuco_clear() |
---|
486 | |
---|
487 | ! Deallocate and reset variables in chemistry module |
---|
488 | CALL chemistry_clear |
---|
489 | |
---|
490 | END SUBROUTINE diffuco_clear |
---|
491 | |
---|
492 | |
---|
493 | !! ================================================================================================================================ |
---|
494 | !! SUBROUTINE : diffuco_aero |
---|
495 | !! |
---|
496 | !>\BRIEF This module first calculates the surface drag |
---|
497 | !! coefficient, for cases in which the surface drag coefficient is NOT provided by the coupled |
---|
498 | !! atmospheric model LMDZ or when the flag ldq_cdrag_from_gcm is set to FALSE |
---|
499 | !! |
---|
500 | !! DESCRIPTION : Computes the surface drag coefficient, for cases |
---|
501 | !! in which it is NOT provided by the coupled atmospheric model LMDZ. The module first uses the |
---|
502 | !! meteorolgical input to calculate the Richardson Number, which is an indicator of atmospheric |
---|
503 | !! stability in the surface layer. The formulation used to find this surface drag coefficient is |
---|
504 | !! dependent on the stability determined. \n |
---|
505 | !! |
---|
506 | !! Designation of wind speed |
---|
507 | !! \latexonly |
---|
508 | !! \input{diffucoaero1.tex} |
---|
509 | !! \endlatexonly |
---|
510 | !! |
---|
511 | !! Calculation of geopotential. This is the definition of Geopotential height (e.g. Jacobson |
---|
512 | !! eqn.4.47, 2005). (required for calculation of the Richardson Number) |
---|
513 | !! \latexonly |
---|
514 | !! \input{diffucoaero2.tex} |
---|
515 | !! \endlatexonly |
---|
516 | !! |
---|
517 | !! \latexonly |
---|
518 | !! \input{diffucoaero3.tex} |
---|
519 | !! \endlatexonly |
---|
520 | !! |
---|
521 | !! Calculation of the virtual air temperature at the surface (required for calculation |
---|
522 | !! of the Richardson Number) |
---|
523 | !! \latexonly |
---|
524 | !! \input{diffucoaero4.tex} |
---|
525 | !! \endlatexonly |
---|
526 | !! |
---|
527 | !! Calculation of the virtual surface temperature (required for calculation of th |
---|
528 | !! Richardson Number) |
---|
529 | !! \latexonly |
---|
530 | !! \input{diffucoaero5.tex} |
---|
531 | !! \endlatexonly |
---|
532 | !! |
---|
533 | !! Calculation of the squared wind shear (required for calculation of the Richardson |
---|
534 | !! Number) |
---|
535 | !! \latexonly |
---|
536 | !! \input{diffucoaero6.tex} |
---|
537 | !! \endlatexonly |
---|
538 | !! |
---|
539 | !! Calculation of the Richardson Number. The Richardson Number is defined as the ratio |
---|
540 | !! of potential to kinetic energy, or, in the context of atmospheric science, of the |
---|
541 | !! generation of energy by wind shear against consumption |
---|
542 | !! by static stability and is an indicator of flow stability (i.e. for when laminar flow |
---|
543 | !! becomes turbulent and vise versa). It is approximated using the expression below: |
---|
544 | !! \latexonly |
---|
545 | !! \input{diffucoaero7.tex} |
---|
546 | !! \endlatexonly |
---|
547 | !! |
---|
548 | !! The Richardson Number hence calculated is subject to a minimum value: |
---|
549 | !! \latexonly |
---|
550 | !! \input{diffucoaero8.tex} |
---|
551 | !! \endlatexonly |
---|
552 | !! |
---|
553 | !! Computing the drag coefficient. We add the add the height of the vegetation to the |
---|
554 | !! level height to take into account that the level 'seen' by the vegetation is actually |
---|
555 | !! the top of the vegetation. Then we we can subtract the displacement height. |
---|
556 | !! \latexonly |
---|
557 | !! \input{diffucoaero9.tex} |
---|
558 | !! \endlatexonly |
---|
559 | !! |
---|
560 | !! For the stable case (i.e $R_i$ $\geq$ 0) |
---|
561 | !! \latexonly |
---|
562 | !! \input{diffucoaero10.tex} |
---|
563 | !! \endlatexonly |
---|
564 | !! |
---|
565 | !! \latexonly |
---|
566 | !! \input{diffucoaero11.tex} |
---|
567 | !! \endlatexonly |
---|
568 | !! |
---|
569 | !! For the unstable case (i.e. $R_i$ < 0) |
---|
570 | !! \latexonly |
---|
571 | !! \input{diffucoaero12.tex} |
---|
572 | !! \endlatexonly |
---|
573 | !! |
---|
574 | !! \latexonly |
---|
575 | !! \input{diffucoaero13.tex} |
---|
576 | !! \endlatexonly |
---|
577 | !! |
---|
578 | !! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere. |
---|
579 | !! To prevent this, a minimum limit to the drag coefficient is defined as: |
---|
580 | !! |
---|
581 | !! \latexonly |
---|
582 | !! \input{diffucoaero14.tex} |
---|
583 | !! \endlatexonly |
---|
584 | !! |
---|
585 | !! RECENT CHANGE(S): None |
---|
586 | !! |
---|
587 | !! MAIN OUTPUT VARIABLE(S): q_cdrag |
---|
588 | !! |
---|
589 | !! REFERENCE(S) : |
---|
590 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
591 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
592 | !! Circulation Model. Journal of Climate, 6, pp.248-273 |
---|
593 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation |
---|
594 | !! sur le cycle de l'eau, PhD Thesis, available from: |
---|
595 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf |
---|
596 | !! - Jacobson M.Z., Fundamentals of Atmospheric Modeling (2nd Edition), published Cambridge |
---|
597 | !! University Press, ISBN 0-521-54865-9 |
---|
598 | !! |
---|
599 | !! FLOWCHART : |
---|
600 | !! \latexonly |
---|
601 | !! \includegraphics[scale=0.5]{diffuco_aero_flowchart.png} |
---|
602 | !! \endlatexonly |
---|
603 | !! \n |
---|
604 | !_ ================================================================================================================================ |
---|
605 | |
---|
606 | SUBROUTINE diffuco_aero (kjpindex, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, & |
---|
607 | qsurf, qair, snow, q_cdrag) |
---|
608 | |
---|
609 | !! 0. Variable and parameter declaration |
---|
610 | |
---|
611 | !! 0.1 Input variables |
---|
612 | |
---|
613 | INTEGER(i_std), INTENT(in) :: kjpindex, kjit !! Domain size |
---|
614 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1}) |
---|
615 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1}) |
---|
616 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first atmospheric layer (m) |
---|
617 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0h !! Surface roughness Length for heat (m) |
---|
618 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0m !! Surface roughness Length for momentum (m) |
---|
619 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective roughness height (m) |
---|
620 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Ground temperature (K) |
---|
621 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature (K) |
---|
622 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! near surface specific air humidity (kg kg^{-1}) |
---|
623 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1}) |
---|
624 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg) |
---|
625 | |
---|
626 | !! 0.2 Output variables |
---|
627 | |
---|
628 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: q_cdrag !! Surface drag coefficient (-) |
---|
629 | |
---|
630 | !! 0.3 Modified variables |
---|
631 | |
---|
632 | !! 0.4 Local variables |
---|
633 | |
---|
634 | INTEGER(i_std) :: ji, jv |
---|
635 | REAL(r_std) :: speed, zg, zdphi, ztvd, ztvs, zdu2 |
---|
636 | REAL(r_std) :: zri, cd_neut, zscf, cd_tmp |
---|
637 | !_ ================================================================================================================================ |
---|
638 | |
---|
639 | !! 1. Initialisation |
---|
640 | |
---|
641 | ! test if we have to work with q_cdrag or to calcul it |
---|
642 | DO ji=1,kjpindex |
---|
643 | |
---|
644 | !! 1a).1 Designation of wind speed |
---|
645 | !! \latexonly |
---|
646 | !! \input{diffucoaero1.tex} |
---|
647 | !! \endlatexonly |
---|
648 | speed = wind(ji) |
---|
649 | |
---|
650 | !! 1a).2 Calculation of geopotentiel |
---|
651 | !! This is the definition of Geopotential height (e.g. Jacobson eqn.4.47, 2005). (required |
---|
652 | !! for calculation of the Richardson Number) |
---|
653 | !! \latexonly |
---|
654 | !! \input{diffucoaero2.tex} |
---|
655 | !! \endlatexonly |
---|
656 | zg = zlev(ji) * cte_grav |
---|
657 | |
---|
658 | !! \latexonly |
---|
659 | !! \input{diffucoaero3.tex} |
---|
660 | !! \endlatexonly |
---|
661 | zdphi = zg/cp_air |
---|
662 | |
---|
663 | !! 1a).3 Calculation of the virtual air temperature at the surface |
---|
664 | !! required for calculation of the Richardson Number |
---|
665 | !! \latexonly |
---|
666 | !! \input{diffucoaero4.tex} |
---|
667 | !! \endlatexonly |
---|
668 | ztvd = (temp_air(ji) + zdphi / (un + rvtmp2 * qair(ji))) * (un + retv * qair(ji)) |
---|
669 | |
---|
670 | !! 1a).4 Calculation of the virtual surface temperature |
---|
671 | !! required for calculation of the Richardson Number |
---|
672 | !! \latexonly |
---|
673 | !! \input{diffucoaero5.tex} |
---|
674 | !! \endlatexonly |
---|
675 | ztvs = temp_sol(ji) * (un + retv * qsurf(ji)) |
---|
676 | |
---|
677 | !! 1a).5 Calculation of the squared wind shear |
---|
678 | !! required for calculation of the Richardson Number |
---|
679 | !! \latexonly |
---|
680 | !! \input{diffucoaero6.tex} |
---|
681 | !! \endlatexonly |
---|
682 | zdu2 = MAX(cepdu2,speed**2) |
---|
683 | |
---|
684 | !! 1a).6 Calculation of the Richardson Number |
---|
685 | !! The Richardson Number is defined as the ratio of potential to kinetic energy, or, in the |
---|
686 | !! context of atmospheric science, of the generation of energy by wind shear against consumption |
---|
687 | !! by static stability and is an indicator of flow stability (i.e. for when laminar flow |
---|
688 | !! becomes turbulent and vise versa).\n |
---|
689 | !! It is approximated using the expression below: |
---|
690 | !! \latexonly |
---|
691 | !! \input{diffucoaero7.tex} |
---|
692 | !! \endlatexonly |
---|
693 | zri = zg * (ztvd - ztvs) / (zdu2 * ztvd) |
---|
694 | |
---|
695 | !! The Richardson Number hence calculated is subject to a minimum value: |
---|
696 | !! \latexonly |
---|
697 | !! \input{diffucoaero8.tex} |
---|
698 | !! \endlatexonly |
---|
699 | zri = MAX(MIN(zri,5.),-5.) |
---|
700 | |
---|
701 | !! 1a).7 Computing the drag coefficient |
---|
702 | !! We add the add the height of the vegetation to the level height to take into account |
---|
703 | !! that the level 'seen' by the vegetation is actually the top of the vegetation. Then we |
---|
704 | !! we can subtract the displacement height. |
---|
705 | !! \latexonly |
---|
706 | !! \input{diffucoaero9.tex} |
---|
707 | !! \endlatexonly |
---|
708 | |
---|
709 | !! 7.0 Snow smoothering |
---|
710 | !! Snow induces low levels of turbulence. |
---|
711 | !! Sensible heat fluxes can therefore be reduced of ~1/3. Pomeroy et al., 1998 |
---|
712 | cd_neut = ct_karman ** 2. / ( LOG( (zlev(ji) + roughheight(ji)) / z0m(ji) ) * LOG( (zlev(ji) + roughheight(ji)) / z0h(ji) ) ) |
---|
713 | |
---|
714 | !! 1a).7.1 - for the stable case (i.e $R_i$ $\geq$ 0) |
---|
715 | IF (zri .GE. zero) THEN |
---|
716 | |
---|
717 | !! \latexonly |
---|
718 | !! \input{diffucoaero10.tex} |
---|
719 | !! \endlatexonly |
---|
720 | zscf = SQRT(un + cd * ABS(zri)) |
---|
721 | |
---|
722 | !! \latexonly |
---|
723 | !! \input{diffucoaero11.tex} |
---|
724 | !! \endlatexonly |
---|
725 | cd_tmp=cd_neut/(un + trois * cb * zri * zscf) |
---|
726 | ELSE |
---|
727 | |
---|
728 | !! 1a).7.2 - for the unstable case (i.e. $R_i$ < 0) |
---|
729 | !! \latexonly |
---|
730 | !! \input{diffucoaero12.tex} |
---|
731 | !! \endlatexonly |
---|
732 | zscf = un / (un + trois * cb * cc * cd_neut * SQRT(ABS(zri) * & |
---|
733 | & ((zlev(ji) + roughheight(ji)) / z0m(ji)))) |
---|
734 | |
---|
735 | !! \latexonly |
---|
736 | !! \input{diffucoaero13.tex} |
---|
737 | !! \endlatexonly |
---|
738 | cd_tmp=cd_neut * (un - trois * cb * zri * zscf) |
---|
739 | ENDIF |
---|
740 | |
---|
741 | !! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere. |
---|
742 | !! To prevent this, a minimum limit to the drag coefficient is defined as: |
---|
743 | |
---|
744 | !! \latexonly |
---|
745 | !! \input{diffucoaero14.tex} |
---|
746 | !! \endlatexonly |
---|
747 | !! |
---|
748 | q_cdrag(ji) = MAX(cd_tmp, min_qc/MAX(speed,min_wind)) |
---|
749 | |
---|
750 | ! In some situations it might be useful to give an upper limit on the cdrag as well. |
---|
751 | ! The line here should then be uncommented. |
---|
752 | !q_cdrag(ji) = MIN(q_cdrag(ji), 0.5/MAX(speed,min_wind)) |
---|
753 | |
---|
754 | END DO |
---|
755 | |
---|
756 | IF (printlev>=3) WRITE (numout,*) ' not ldqcdrag_from_gcm : diffuco_aero done ' |
---|
757 | |
---|
758 | END SUBROUTINE diffuco_aero |
---|
759 | |
---|
760 | |
---|
761 | !! ================================================================================================================================ |
---|
762 | !! SUBROUTINE : diffuco_snow |
---|
763 | !! |
---|
764 | !>\BRIEF This subroutine computes the beta coefficient for snow sublimation. |
---|
765 | !! |
---|
766 | !! DESCRIPTION : This routine computes beta coefficient for snow sublimation, which |
---|
767 | !! integrates the snow on both vegetation and other surface types (e.g. ice, lakes, |
---|
768 | !! cities etc.) \n |
---|
769 | !! |
---|
770 | !! A critical depth of snow (snowcri) is defined to calculate the fraction of each grid-cell |
---|
771 | !! that is covered with snow (snow/snowcri) while the remaining part is snow-free. |
---|
772 | !! We also carry out a first calculation of sublimation (subtest) to lower down the beta |
---|
773 | !! coefficient if necessary (if subtest > snow). This is a predictor-corrector test. |
---|
774 | !! |
---|
775 | !! RECENT CHANGE(S): None |
---|
776 | !! |
---|
777 | !! MAIN OUTPUT VARIABLE(S): ::vbeta1 |
---|
778 | !! |
---|
779 | !! REFERENCE(S) : |
---|
780 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
781 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
782 | !! Circulation Model. Journal of Climate, 6, pp. 248-273 |
---|
783 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation |
---|
784 | !! sur le cycle de l'eau, PhD Thesis, available from: |
---|
785 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf |
---|
786 | !! |
---|
787 | !! FLOWCHART : None |
---|
788 | !! \n |
---|
789 | !_ ================================================================================================================================ |
---|
790 | |
---|
791 | SUBROUTINE diffuco_snow (kjpindex, qair, qsatt, rau, u, v,q_cdrag, & |
---|
792 | & snow, frac_nobio, totfrac_nobio, snow_nobio, frac_snow_veg, frac_snow_nobio, & |
---|
793 | vbeta1) |
---|
794 | |
---|
795 | !! 0. Variable and parameter declaration |
---|
796 | |
---|
797 | !! 0.1 Input variables |
---|
798 | |
---|
799 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
800 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1}) |
---|
801 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1}) |
---|
802 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3}) |
---|
803 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1}) |
---|
804 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1}) |
---|
805 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
806 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg m^{-2}) |
---|
807 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, cities etc. (-) |
---|
808 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice, lakes, cities etc. (-) |
---|
809 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice, lakes, cities etc. (-) |
---|
810 | REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area |
---|
811 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio!! Snow cover fraction on non-vegeted area |
---|
812 | |
---|
813 | !! 0.2 Output variables |
---|
814 | |
---|
815 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation (dimensionless ratio) |
---|
816 | |
---|
817 | !! 0.3 Modified variables |
---|
818 | |
---|
819 | !! 0.4 Local variables |
---|
820 | |
---|
821 | REAL(r_std) :: subtest !! Sublimation for test (kg m^{-2}) |
---|
822 | REAL(r_std) :: zrapp !! Modified factor (ratio) |
---|
823 | REAL(r_std) :: speed !! Wind speed (m s^{-1}) |
---|
824 | REAL(r_std) :: vbeta1_add !! Beta for sublimation (ratio) |
---|
825 | INTEGER(i_std) :: ji, jv !! Indices (-) |
---|
826 | !_ ================================================================================================================================ |
---|
827 | |
---|
828 | !! 1. Calculate beta coefficient for snow sublimation on the vegetation\n |
---|
829 | |
---|
830 | DO ji=1,kjpindex ! Loop over # pixels - domain size |
---|
831 | |
---|
832 | ! Fraction of mesh that can sublimate snow |
---|
833 | vbeta1(ji) = (un - totfrac_nobio(ji)) * frac_snow_veg(ji) |
---|
834 | |
---|
835 | ! Limitation of sublimation in case of snow amounts smaller than the atmospheric demand. |
---|
836 | speed = MAX(min_wind, wind(ji)) |
---|
837 | |
---|
838 | subtest = dt_sechiba * vbeta1(ji) * speed * q_cdrag(ji) * rau(ji) * & |
---|
839 | & ( qsatt(ji) - qair(ji) ) |
---|
840 | |
---|
841 | IF ( subtest .GT. min_sechiba ) THEN |
---|
842 | zrapp = snow(ji) / subtest |
---|
843 | IF ( zrapp .LT. un ) THEN |
---|
844 | vbeta1(ji) = vbeta1(ji) * zrapp |
---|
845 | ENDIF |
---|
846 | ENDIF |
---|
847 | |
---|
848 | END DO ! Loop over # pixels - domain size |
---|
849 | |
---|
850 | !! 2. Add the beta coefficients calculated from other surfaces types (snow on ice,lakes, cities...) |
---|
851 | |
---|
852 | DO jv = 1, nnobio ! Loop over # other surface types |
---|
853 | !!$ ! |
---|
854 | !!$ IF ( jv .EQ. iice ) THEN |
---|
855 | !!$ ! |
---|
856 | !!$ ! Land ice is of course a particular case |
---|
857 | !!$ ! |
---|
858 | !!$ DO ji=1,kjpindex |
---|
859 | !!$ vbeta1(ji) = vbeta1(ji) + frac_nobio(ji,jv) |
---|
860 | !!$ ENDDO |
---|
861 | !!$ ! |
---|
862 | !!$ ELSE |
---|
863 | ! |
---|
864 | DO ji=1,kjpindex ! Loop over # pixels - domain size |
---|
865 | |
---|
866 | vbeta1_add = frac_nobio(ji,jv) * frac_snow_nobio(ji, jv) |
---|
867 | |
---|
868 | ! Limitation of sublimation in case of snow amounts smaller than |
---|
869 | ! the atmospheric demand. |
---|
870 | speed = MAX(min_wind, wind(ji)) |
---|
871 | |
---|
872 | !! Limitation of sublimation by the snow accumulated on the ground |
---|
873 | !! A first approximation is obtained with the old values of |
---|
874 | !! qair and qsol_sat: function of temp-sol and pb. (see call of qsatcalc) |
---|
875 | subtest = dt_sechiba * vbeta1_add * speed * q_cdrag(ji) * rau(ji) * & |
---|
876 | & ( qsatt(ji) - qair(ji) ) |
---|
877 | |
---|
878 | IF ( subtest .GT. min_sechiba ) THEN |
---|
879 | zrapp = snow_nobio(ji,jv) / subtest |
---|
880 | IF ( zrapp .LT. un ) THEN |
---|
881 | vbeta1_add = vbeta1_add * zrapp |
---|
882 | ENDIF |
---|
883 | ENDIF |
---|
884 | |
---|
885 | vbeta1(ji) = vbeta1(ji) + vbeta1_add |
---|
886 | |
---|
887 | ENDDO ! Loop over # pixels - domain size |
---|
888 | |
---|
889 | !!$ ENDIF |
---|
890 | |
---|
891 | ENDDO ! Loop over # other surface types |
---|
892 | |
---|
893 | IF (printlev>=3) WRITE (numout,*) ' diffuco_snow done ' |
---|
894 | |
---|
895 | END SUBROUTINE diffuco_snow |
---|
896 | |
---|
897 | |
---|
898 | !! ================================================================================================================================ |
---|
899 | !! SUBROUTINE : diffuco_flood |
---|
900 | !! |
---|
901 | !>\BRIEF This routine computes partial beta coefficient : floodplains |
---|
902 | !! |
---|
903 | !! DESCRIPTION : |
---|
904 | !! |
---|
905 | !! RECENT CHANGE(S): None |
---|
906 | !! |
---|
907 | !! MAIN OUTPUT VARIABLE(S) : vbeta5 |
---|
908 | !! |
---|
909 | !! REFERENCE(S) : None |
---|
910 | !! |
---|
911 | !! FLOWCHART : None |
---|
912 | !! \n |
---|
913 | !_ ================================================================================================================================ |
---|
914 | |
---|
915 | SUBROUTINE diffuco_flood (kjpindex, qair, qsatt, rau, u, v, q_cdrag, evapot, evapot_corr, & |
---|
916 | & flood_frac, flood_res, vbeta5) |
---|
917 | |
---|
918 | ! interface description |
---|
919 | ! input scalar |
---|
920 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size |
---|
921 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity |
---|
922 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity |
---|
923 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density |
---|
924 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed |
---|
925 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed |
---|
926 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
927 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_res !! water mass in flood reservoir |
---|
928 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_frac !! fraction of floodplains |
---|
929 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Potential evaporation |
---|
930 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr!! Potential evaporation2 |
---|
931 | ! output fields |
---|
932 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta5 !! Beta for floodplains |
---|
933 | |
---|
934 | ! local declaration |
---|
935 | REAL(r_std) :: subtest, zrapp, speed |
---|
936 | INTEGER(i_std) :: ji, jv |
---|
937 | |
---|
938 | !_ ================================================================================================================================ |
---|
939 | ! |
---|
940 | ! beta coefficient for sublimation for floodplains |
---|
941 | ! |
---|
942 | DO ji=1,kjpindex |
---|
943 | ! |
---|
944 | IF (evapot(ji) .GT. min_sechiba) THEN |
---|
945 | vbeta5(ji) = flood_frac(ji) *evapot_corr(ji)/evapot(ji) |
---|
946 | ELSE |
---|
947 | vbeta5(ji) = flood_frac(ji) |
---|
948 | ENDIF |
---|
949 | ! |
---|
950 | ! -- Limitation of evaporation in case of water amounts smaller than |
---|
951 | ! the atmospheric demand. |
---|
952 | |
---|
953 | ! |
---|
954 | speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) |
---|
955 | ! |
---|
956 | subtest = dt_sechiba * vbeta5(ji) * speed * q_cdrag(ji) * rau(ji) * & |
---|
957 | & ( qsatt(ji) - qair(ji) ) |
---|
958 | ! |
---|
959 | IF ( subtest .GT. min_sechiba ) THEN |
---|
960 | zrapp = flood_res(ji) / subtest |
---|
961 | IF ( zrapp .LT. un ) THEN |
---|
962 | vbeta5(ji) = vbeta5(ji) * zrapp |
---|
963 | ENDIF |
---|
964 | ENDIF |
---|
965 | ! |
---|
966 | END DO |
---|
967 | |
---|
968 | IF (printlev>=3) WRITE (numout,*) ' diffuco_flood done ' |
---|
969 | |
---|
970 | END SUBROUTINE diffuco_flood |
---|
971 | |
---|
972 | |
---|
973 | !! ================================================================================================================================ |
---|
974 | !! SUBROUTINE : diffuco_inter |
---|
975 | !! |
---|
976 | !>\BRIEF This routine computes the partial beta coefficient |
---|
977 | !! for the interception for each type of vegetation |
---|
978 | !! |
---|
979 | !! DESCRIPTION : We first calculate the dry and wet parts of each PFT (wet part = qsintveg/qsintmax). |
---|
980 | !! The former is submitted to transpiration only (vbeta3 coefficient, calculated in |
---|
981 | !! diffuco_trans_co2), while the latter is first submitted to interception loss |
---|
982 | !! (vbeta2 coefficient) and then to transpiration once all the intercepted water has been evaporated |
---|
983 | !! (vbeta23 coefficient). Interception loss is also submitted to a predictor-corrector test, |
---|
984 | !! as for snow sublimation. \n |
---|
985 | !! |
---|
986 | !! \latexonly |
---|
987 | !! \input{diffucointer1.tex} |
---|
988 | !! \endlatexonly |
---|
989 | !! Calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water |
---|
990 | !! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss. |
---|
991 | !! |
---|
992 | !! \latexonly |
---|
993 | !! \input{diffucointer2.tex} |
---|
994 | !! \endlatexonly |
---|
995 | !! |
---|
996 | !! Calculation of $\beta_3$, the canopy transpiration resistance |
---|
997 | !! \latexonly |
---|
998 | !! \input{diffucointer3.tex} |
---|
999 | !! \endlatexonly |
---|
1000 | !! |
---|
1001 | !! We here determine the limitation of interception loss by the water stored on the leaf. |
---|
1002 | !! A first approximation of interception loss is obtained using the old values of |
---|
1003 | !! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc') |
---|
1004 | !! \latexonly |
---|
1005 | !! \input{diffucointer4.tex} |
---|
1006 | !! \endlatexonly |
---|
1007 | !! |
---|
1008 | !! \latexonly |
---|
1009 | !! \input{diffucointer5.tex} |
---|
1010 | !! \endlatexonly |
---|
1011 | !! |
---|
1012 | !! \latexonly |
---|
1013 | !! \input{diffucointer6.tex} |
---|
1014 | !! \endlatexonly |
---|
1015 | !! |
---|
1016 | !! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction |
---|
1017 | !! 'zqsvegrap'. |
---|
1018 | !! \latexonly |
---|
1019 | !! \input{diffucointer7.tex} |
---|
1020 | !! \endlatexonly |
---|
1021 | !! |
---|
1022 | !! RECENT CHANGE(S): None |
---|
1023 | !! |
---|
1024 | !! MAIN OUTPUT VARIABLE(S): ::vbeta2, ::vbeta23 |
---|
1025 | !! |
---|
1026 | !! REFERENCE(S) : |
---|
1027 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
1028 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
1029 | !! Circulation Model. Journal of Climate, 6, pp. 248-273 |
---|
1030 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation |
---|
1031 | !! sur le cycle de l'eau, PhD Thesis, available from: |
---|
1032 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf |
---|
1033 | !! - Perrier, A, 1975. Etude physique de l'évaporation dans les conditions naturelles. Annales |
---|
1034 | !! Agronomiques, 26(1-18): pp. 105-123, pp. 229-243 |
---|
1035 | !! |
---|
1036 | !! FLOWCHART : None |
---|
1037 | !! \n |
---|
1038 | !_ ================================================================================================================================ |
---|
1039 | |
---|
1040 | SUBROUTINE diffuco_inter (kjpindex, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, & |
---|
1041 | & qsintveg, qsintmax, rstruct, vbeta2, vbeta23) |
---|
1042 | |
---|
1043 | !! 0 Variable and parameter declaration |
---|
1044 | |
---|
1045 | !! 0.1 Input variables |
---|
1046 | |
---|
1047 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
1048 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1}) |
---|
1049 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1}) |
---|
1050 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3}) |
---|
1051 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1}) |
---|
1052 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1}) |
---|
1053 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
1054 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1) |
---|
1055 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! vegetation fraction for each type (fraction) |
---|
1056 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2}) |
---|
1057 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2}) |
---|
1058 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! architectural resistance (s m^{-1}) |
---|
1059 | |
---|
1060 | !! 0.2 Output variables |
---|
1061 | |
---|
1062 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss (-) |
---|
1063 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta23 !! Beta for fraction of wetted foliage that will |
---|
1064 | !! transpire (-) |
---|
1065 | |
---|
1066 | !! 0.4 Local variables |
---|
1067 | |
---|
1068 | INTEGER(i_std) :: ji, jv !! (-), (-) |
---|
1069 | REAL(r_std) :: zqsvegrap, ziltest, zrapp, speed !! |
---|
1070 | !_ ================================================================================================================================ |
---|
1071 | |
---|
1072 | !! 1. Initialize |
---|
1073 | |
---|
1074 | vbeta2(:,:) = zero |
---|
1075 | vbeta23(:,:) = zero |
---|
1076 | |
---|
1077 | !! 2. The beta coefficient for interception by vegetation. |
---|
1078 | |
---|
1079 | DO jv = 2,nvm |
---|
1080 | |
---|
1081 | DO ji=1,kjpindex |
---|
1082 | |
---|
1083 | IF (veget(ji,jv) .GT. min_sechiba .AND. qsintveg(ji,jv) .GT. zero ) THEN |
---|
1084 | |
---|
1085 | zqsvegrap = zero |
---|
1086 | IF (qsintmax(ji,jv) .GT. min_sechiba ) THEN |
---|
1087 | |
---|
1088 | !! \latexonly |
---|
1089 | !! \input{diffucointer1.tex} |
---|
1090 | !! \endlatexonly |
---|
1091 | !! |
---|
1092 | !! We calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water |
---|
1093 | !! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss. |
---|
1094 | !! |
---|
1095 | zqsvegrap = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv)) |
---|
1096 | END IF |
---|
1097 | |
---|
1098 | !! \latexonly |
---|
1099 | !! \input{diffucointer2.tex} |
---|
1100 | !! \endlatexonly |
---|
1101 | speed = MAX(min_wind, wind(ji)) |
---|
1102 | |
---|
1103 | !! Calculation of $\beta_3$, the canopy transpiration resistance |
---|
1104 | !! \latexonly |
---|
1105 | !! \input{diffucointer3.tex} |
---|
1106 | !! \endlatexonly |
---|
1107 | vbeta2(ji,jv) = veget(ji,jv) * zqsvegrap * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv))) |
---|
1108 | |
---|
1109 | !! We here determine the limitation of interception loss by the water stored on the leaf. |
---|
1110 | !! A first approximation of interception loss is obtained using the old values of |
---|
1111 | !! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc') |
---|
1112 | !! \latexonly |
---|
1113 | !! \input{diffucointer4.tex} |
---|
1114 | !! \endlatexonly |
---|
1115 | ziltest = dt_sechiba * vbeta2(ji,jv) * speed * q_cdrag(ji) * rau(ji) * & |
---|
1116 | & ( qsatt(ji) - qair(ji) ) |
---|
1117 | |
---|
1118 | IF ( ziltest .GT. min_sechiba ) THEN |
---|
1119 | |
---|
1120 | !! \latexonly |
---|
1121 | !! \input{diffucointer5.tex} |
---|
1122 | !! \endlatexonly |
---|
1123 | zrapp = qsintveg(ji,jv) / ziltest |
---|
1124 | IF ( zrapp .LT. un ) THEN |
---|
1125 | |
---|
1126 | !! \latexonly |
---|
1127 | !! \input{diffucointer6.tex} |
---|
1128 | !! \endlatexonly |
---|
1129 | !! |
---|
1130 | !! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction |
---|
1131 | !! 'zqsvegrap'. |
---|
1132 | IF ( humrel(ji,jv) >= min_sechiba ) THEN |
---|
1133 | vbeta23(ji,jv) = MAX(vbeta2(ji,jv) - vbeta2(ji,jv) * zrapp, zero) |
---|
1134 | ELSE |
---|
1135 | ! We don't want transpiration when the soil cannot deliver it |
---|
1136 | vbeta23(ji,jv) = zero |
---|
1137 | ENDIF |
---|
1138 | |
---|
1139 | !! \latexonly |
---|
1140 | !! \input{diffucointer7.tex} |
---|
1141 | !! \endlatexonly |
---|
1142 | vbeta2(ji,jv) = vbeta2(ji,jv) * zrapp |
---|
1143 | ENDIF |
---|
1144 | ENDIF |
---|
1145 | END IF |
---|
1146 | ! ! Autre formulation possible pour l'evaporation permettant une transpiration sur tout le feuillage |
---|
1147 | ! !commenter si formulation Nathalie sinon Tristan |
---|
1148 | ! speed = MAX(min_wind, wind(ji)) |
---|
1149 | ! |
---|
1150 | ! vbeta23(ji,jv) = MAX(zero, veget(ji,jv) * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv))) - vbeta2(ji,jv)) |
---|
1151 | |
---|
1152 | END DO |
---|
1153 | |
---|
1154 | END DO |
---|
1155 | |
---|
1156 | IF (printlev>=3) WRITE (numout,*) ' diffuco_inter done ' |
---|
1157 | |
---|
1158 | END SUBROUTINE diffuco_inter |
---|
1159 | |
---|
1160 | |
---|
1161 | !! ============================================================================================================================== |
---|
1162 | !! SUBROUTINE : diffuco_bare |
---|
1163 | !! |
---|
1164 | !>\BRIEF This routine computes the partial beta coefficient corresponding to |
---|
1165 | !! bare soil |
---|
1166 | !! |
---|
1167 | !! DESCRIPTION : Bare soil evaporation is submitted to a maximum possible flow (evap_bare_lim) |
---|
1168 | !! |
---|
1169 | !! Calculation of wind speed |
---|
1170 | !! \latexonly |
---|
1171 | !! \input{diffucobare1.tex} |
---|
1172 | !! \endlatexonly |
---|
1173 | !! |
---|
1174 | !! The calculation of $\beta_4$ |
---|
1175 | !! \latexonly |
---|
1176 | !! \input{diffucobare2.tex} |
---|
1177 | !! \endlatexonly |
---|
1178 | !! |
---|
1179 | !! RECENT CHANGE(S): None |
---|
1180 | !! |
---|
1181 | !! MAIN OUTPUT VARIABLE(S): ::vbeta4 |
---|
1182 | !! |
---|
1183 | !! REFERENCE(S) : |
---|
1184 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
1185 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
1186 | !! Circulation Model. Journal of Climate, 6, pp.248-273 |
---|
1187 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation |
---|
1188 | !! sur le cycle de l'eau, PhD Thesis, available from: |
---|
1189 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf |
---|
1190 | !! |
---|
1191 | !! FLOWCHART : None |
---|
1192 | !! \n |
---|
1193 | !_ ================================================================================================================================ |
---|
1194 | |
---|
1195 | SUBROUTINE diffuco_bare (kjpindex, tot_bare_soil, veget_max, evap_bare_lim, evap_bare_lim_ns, vbeta2, vbeta3, vbeta4) |
---|
1196 | |
---|
1197 | !! 0. Variable and parameter declaration |
---|
1198 | |
---|
1199 | !! 0.1 Input variables |
---|
1200 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
1201 | REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction |
---|
1202 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty) |
---|
1203 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Beta for Interception |
---|
1204 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Beta for Transpiration |
---|
1205 | |
---|
1206 | !! 0.2 Output variables |
---|
1207 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation (-) |
---|
1208 | |
---|
1209 | !! 0.3 Modified variables |
---|
1210 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evap_bare_lim !! limiting factor for bare soil evaporation |
---|
1211 | !! when the 11-layer hydrology is used (-) |
---|
1212 | REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout) :: evap_bare_lim_ns !! limiting factor for bare soil evaporation |
---|
1213 | !! when the 11-layer hydrology is used (-) |
---|
1214 | !! 0.4 Local variables |
---|
1215 | INTEGER(i_std) :: ji |
---|
1216 | REAL(r_std), DIMENSION(kjpindex) :: vegtot |
---|
1217 | !_ ================================================================================================================================ |
---|
1218 | |
---|
1219 | !! 1. Calculation of the soil resistance and the beta (beta_4) for bare soil |
---|
1220 | |
---|
1221 | ! To use the new version of hydrol_split_soil and ensure the water conservation, we must have throughout sechiba at all time: |
---|
1222 | ! a) evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)*soiltile(ji,:)*vegtot(ji)) |
---|
1223 | ! b) all the terms (vbeta4, evap_bare_lim, evap_bare_lim_ns) =0 if evap_bare_lim(ji) LE min_sechiba |
---|
1224 | ! This must also be kept true in diffuco_comb |
---|
1225 | |
---|
1226 | DO ji = 1, kjpindex |
---|
1227 | |
---|
1228 | ! The limitation by 1-beta2-beta3 is due to the fact that evaporation under vegetation is possible |
---|
1229 | !! \latexonly |
---|
1230 | !! \input{diffucobare3.tex} |
---|
1231 | !! \endlatexonly |
---|
1232 | |
---|
1233 | IF ( (evap_bare_lim(ji) .GT. min_sechiba) .AND. & |
---|
1234 | ! in this case we can't have vegtot LE min_sechina, cf hydrol_soil |
---|
1235 | (un - SUM(vbeta2(ji,:)+vbeta3(ji,:)) .GT. min_sechiba) ) THEN |
---|
1236 | ! eventually, none of the left-hand term is close to zero |
---|
1237 | |
---|
1238 | vegtot(ji) = SUM(veget_max(ji,:)) |
---|
1239 | IF (evap_bare_lim(ji) < (un - SUM(vbeta2(ji,:)+vbeta3(ji,:)))) THEN |
---|
1240 | ! Standard case |
---|
1241 | vbeta4(ji) = evap_bare_lim(ji) |
---|
1242 | ELSE |
---|
1243 | vbeta4(ji) = un - SUM(vbeta2(ji,:)+vbeta3(ji,:)) |
---|
1244 | |
---|
1245 | ! We now have to redefine evap_bare_lim & evap_bare_lim_ns |
---|
1246 | IF (evap_bare_lim(ji) .GT. min_sechiba) THEN |
---|
1247 | evap_bare_lim_ns(ji,:) = evap_bare_lim_ns(ji,:) * vbeta4(ji) / evap_bare_lim(ji) |
---|
1248 | ELSE ! we must re-invent evap_bare_lim_ns => uniform across soiltiles |
---|
1249 | evap_bare_lim_ns(ji,:) = tot_bare_soil(ji)/vegtot(ji) |
---|
1250 | ENDIF |
---|
1251 | |
---|
1252 | evap_bare_lim(ji) = vbeta4(ji) |
---|
1253 | ! consistent with evap_bare_lim(ji) = |
---|
1254 | ! SUM(evap_bare_lim_ns(ji,:)*soiltile(ji,:)*vegtot(ji)) |
---|
1255 | ! as SUM(soiltile(ji,:)) = 1 |
---|
1256 | END IF |
---|
1257 | |
---|
1258 | ELSE ! instead of having very small vbeta4, we set everything to zero |
---|
1259 | vbeta4(ji) = zero |
---|
1260 | evap_bare_lim(ji) = zero |
---|
1261 | evap_bare_lim_ns(ji,:) = zero |
---|
1262 | ENDIF |
---|
1263 | |
---|
1264 | END DO |
---|
1265 | |
---|
1266 | IF (printlev>=3) WRITE (numout,*) ' diffuco_bare done ' |
---|
1267 | |
---|
1268 | END SUBROUTINE diffuco_bare |
---|
1269 | |
---|
1270 | |
---|
1271 | !! ============================================================================================================================== |
---|
1272 | !! SUBROUTINE : diffuco_trans_co2 |
---|
1273 | !! |
---|
1274 | !>\BRIEF This subroutine computes carbon assimilation and stomatal |
---|
1275 | !! conductance, following respectively Farqhuar et al. (1980) and Ball et al. (1987). |
---|
1276 | !! |
---|
1277 | !! DESCRIPTION :\n |
---|
1278 | !! *** General:\n |
---|
1279 | !! The equations are different depending on the photosynthesis mode (C3 versus C4).\n |
---|
1280 | !! Assimilation and conductance are computed over 20 levels of LAI and then |
---|
1281 | !! integrated at the canopy level.\n |
---|
1282 | !! This routine also computes partial beta coefficient: transpiration for each |
---|
1283 | !! type of vegetation.\n |
---|
1284 | !! There is a main loop on the PFTs, then inner loops on the points where |
---|
1285 | !! assimilation has to be calculated.\n |
---|
1286 | !! This subroutine is called at each sechiba time step by sechiba_main.\n |
---|
1287 | !! *** Details: |
---|
1288 | !! - Integration at the canopy level\n |
---|
1289 | !! \latexonly |
---|
1290 | !! \input{diffuco_trans_co2_1.1.tex} |
---|
1291 | !! \endlatexonly |
---|
1292 | !! - Light''s extinction \n |
---|
1293 | !! The available light follows a simple Beer extinction law. |
---|
1294 | !! The extinction coefficients (ext_coef) are PFT-dependant constants and are defined in constant_co2.f90.\n |
---|
1295 | !! \latexonly |
---|
1296 | !! \input{diffuco_trans_co2_1.2.tex} |
---|
1297 | !! \endlatexonly |
---|
1298 | !! - Estimation of relative humidity of air (for calculation of the stomatal conductance)\n |
---|
1299 | !! \latexonly |
---|
1300 | !! \input{diffuco_trans_co2_1.3.tex} |
---|
1301 | !! \endlatexonly |
---|
1302 | !! - Calculation of the water limitation factor\n |
---|
1303 | !! \latexonly |
---|
1304 | !! \input{diffuco_trans_co2_2.1.tex} |
---|
1305 | !! \endlatexonly |
---|
1306 | !! - Calculation of temperature dependent parameters for C4 plants\n |
---|
1307 | !! \latexonly |
---|
1308 | !! \input{diffuco_trans_co2_2.2.tex} |
---|
1309 | !! \endlatexonly |
---|
1310 | !! - Calculation of temperature dependent parameters for C3 plants\n |
---|
1311 | !! \latexonly |
---|
1312 | !! \input{diffuco_trans_co2_2.3.tex} |
---|
1313 | !! \endlatexonly |
---|
1314 | !! - Vmax scaling\n |
---|
1315 | !! Vmax is scaled into the canopy due to reduction of nitrogen |
---|
1316 | !! (Johnson and Thornley,1984).\n |
---|
1317 | !! \latexonly |
---|
1318 | !! \input{diffuco_trans_co2_2.4.1.tex} |
---|
1319 | !! \endlatexonly |
---|
1320 | !! - Assimilation for C4 plants (Collatz et al., 1992)\n |
---|
1321 | !! \latexonly |
---|
1322 | !! \input{diffuco_trans_co2_2.4.2.tex} |
---|
1323 | !! \endlatexonly |
---|
1324 | !! - Assimilation for C3 plants (Farqhuar et al., 1980)\n |
---|
1325 | !! \latexonly |
---|
1326 | !! \input{diffuco_trans_co2_2.4.3.tex} |
---|
1327 | !! \endlatexonly |
---|
1328 | !! - Estimation of the stomatal conductance (Ball et al., 1987)\n |
---|
1329 | !! \latexonly |
---|
1330 | !! \input{diffuco_trans_co2_2.4.4.tex} |
---|
1331 | !! \endlatexonly |
---|
1332 | !! |
---|
1333 | !! RECENT CHANGE(S): |
---|
1334 | !! |
---|
1335 | !! 2018/12 added a nitrogen use efficiency constraint on maximum rate of carboxylation |
---|
1336 | !! |
---|
1337 | !! 2018/11 replaced t2m and q2m by temp_air and qair; prognostic variables on first atmospheric model layer |
---|
1338 | !! |
---|
1339 | !! 2006/06 N. de Noblet |
---|
1340 | !! - addition of q2m and t2m as input parameters for the |
---|
1341 | !! calculation of Rveget |
---|
1342 | !! - introduction of vbeta23 |
---|
1343 | !! |
---|
1344 | !! MAIN OUTPUT VARIABLE(S): beta coefficients, resistances, CO2 intercellular |
---|
1345 | !! concentration |
---|
1346 | !! |
---|
1347 | !! REFERENCE(S) : |
---|
1348 | !! - Ball, J., T. Woodrow, and J. Berry (1987), A model predicting stomatal |
---|
1349 | !! conductance and its contribution to the control of photosynthesis under |
---|
1350 | !! different environmental conditions, Prog. Photosynthesis, 4, 221â 224. |
---|
1351 | !! - Collatz, G., M. Ribas-Carbo, and J. Berry (1992), Coupled photosynthesis |
---|
1352 | !! stomatal conductance model for leaves of C4 plants, Aust. J. Plant Physiol., |
---|
1353 | !! 19, 519â538. |
---|
1354 | !! - Farquhar, G., S. von Caemmener, and J. Berry (1980), A biochemical model of |
---|
1355 | !! photosynthesis CO2 fixation in leaves of C3 species, Planta, 149, 78â90. |
---|
1356 | !! - Johnson, I. R., and J. Thornley (1984), A model of instantaneous and daily |
---|
1357 | !! canopy photosynthesis, J Theor. Biol., 107, 531 545 |
---|
1358 | !! - McMurtrie, R.E., Rook, D.A. and Kelliher, F.M., 1990. Modelling the yield of Pinus radiata on a |
---|
1359 | !! site limited by water and nitrogen. For. Ecol. Manage., 30: 381-413 |
---|
1360 | !! - Bounoua, L., Hall, F. G., Sellers, P. J., Kumar, A., Collatz, G. J., Tucker, C. J., and Imhoff, M. L. (2010), Quantifying the |
---|
1361 | !! negative feedback of vegetation to greenhouse warming: A modeling approach, Geophysical Research Letters, 37, Artn L23701, |
---|
1362 | !! Doi 10.1029/2010gl045338 |
---|
1363 | !! - Bounoua, L., Collatz, G. J., Sellers, P. J., Randall, D. A., Dazlich, D. A., Los, S. O., Berry, J. A., Fung, I., |
---|
1364 | !! Tucker, C. J., Field, C. B., and Jensen, T. G. (1999), Interactions between vegetation and climate: Radiative and physiological |
---|
1365 | !! effects of doubled atmospheric co2, Journal of Climate, 12, 309-324, Doi 10.1175/1520-0442(1999)012<0309:Ibvacr>2.0.Co;2 |
---|
1366 | !! - Sellers, P. J., Bounoua, L., Collatz, G. J., Randall, D. A., Dazlich, D. A., Los, S. O., Berry, J. A., Fung, I., |
---|
1367 | !! Tucker, C. J., Field, C. B., and Jensen, T. G. (1996), Comparison of radiative and physiological effects of doubled atmospheric |
---|
1368 | !! co2 on climate, Science, 271, 1402-1406, DOI 10.1126/science.271.5254.1402 |
---|
1369 | !! - Lewis, J. D., Ward, J. K., and Tissue, D. T. (2010), Phosphorus supply drives nonlinear responses of cottonwood |
---|
1370 | !! (populus deltoides) to increases in co2 concentration from glacial to future concentrations, New Phytologist, 187, 438-448, |
---|
1371 | !! DOI 10.1111/j.1469-8137.2010.03307.x |
---|
1372 | !! - Kattge, J., Knorr, W., Raddatz, T., and Wirth, C. (2009), Quantifying photosynthetic capacity and its relationship to leaf |
---|
1373 | !! nitrogen content for global-scale terrestrial biosphere models, Global Change Biology, 15, 976-991, |
---|
1374 | !! DOI 10.1111/j.1365-2486.2008.01744.x |
---|
1375 | !! |
---|
1376 | !! FLOWCHART : None |
---|
1377 | !! \n |
---|
1378 | !_ ================================================================================================================================ |
---|
1379 | |
---|
1380 | SUBROUTINE diffuco_trans_co2 (kjpindex, swdown, pb, qsurf, qair, temp_air, temp_growth, rau, u, v, q_cdrag, humrel, & |
---|
1381 | assim_param, Ca, & |
---|
1382 | veget, veget_max, lai, qsintveg, qsintmax, vbeta3, vbeta3pot, rveget, rstruct, & |
---|
1383 | cimean, gsmean, gpp, & |
---|
1384 | vbeta23, hist_id, indexveg, indexlai, index, kjit, cim) |
---|
1385 | |
---|
1386 | ! |
---|
1387 | !! 0. Variable and parameter declaration |
---|
1388 | ! |
---|
1389 | |
---|
1390 | ! |
---|
1391 | !! 0.1 Input variables |
---|
1392 | ! |
---|
1393 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
1394 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Downwelling short wave flux |
---|
1395 | !! @tex ($W m^{-2}$) @endtex |
---|
1396 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa) |
---|
1397 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Near surface specific humidity |
---|
1398 | !! @tex ($kg kg^{-1}$) @endtex |
---|
1399 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Specific humidity at first atmospheric model layer |
---|
1400 | !! @tex ($kg kg^{-1}$) @endtex |
---|
1401 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature at first atmospheric model layer (K) |
---|
1402 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_growth !! Growth temperature (°C) - Is equal to t2m_month |
---|
1403 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! air density @tex ($kg m^{-3}$) @endtex |
---|
1404 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed |
---|
1405 | !! @tex ($m s^{-1}$) @endtex |
---|
1406 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed |
---|
1407 | !! @tex ($m s^{-1}$) @endtex |
---|
1408 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
1409 | REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param !! referenced vcmax, nue, leaf N |
---|
1410 | !! for photosynthesis |
---|
1411 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: Ca !! CO2 concentration inside the canopy |
---|
1412 | !! @tex ($\mu mol mol^{-1}$) @endtex |
---|
1413 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (0-1,unitless) |
---|
1414 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Coverage fraction of vegetation for each PFT |
---|
1415 | !! depending on LAI (0-1, unitless) |
---|
1416 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Maximum vegetation fraction of each PFT inside |
---|
1417 | !! the grid box (0-1, unitless) |
---|
1418 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index @tex ($m^2 m^{-2}$) @endtex |
---|
1419 | !! @tex ($m s^{-1}$) @endtex |
---|
1420 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception |
---|
1421 | !! @tex ($kg m^{-2}$) @endte |
---|
1422 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation |
---|
1423 | !! @tex ($kg m^{-2}$) @endtex |
---|
1424 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23 !! Beta for fraction of wetted foliage that will |
---|
1425 | !! transpire (unitless) |
---|
1426 | INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (-) |
---|
1427 | INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map (-) |
---|
1428 | INTEGER(i_std),DIMENSION (kjpindex*(nlai+1)), INTENT (in) :: indexlai !! Indeces of the points on the 3D map |
---|
1429 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-) |
---|
1430 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-) |
---|
1431 | ! |
---|
1432 | !! 0.2 Output variables |
---|
1433 | ! |
---|
1434 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for Transpiration (unitless) |
---|
1435 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3pot !! Beta for Potential Transpiration |
---|
1436 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! stomatal resistance of vegetation |
---|
1437 | !! @tex ($s m^{-1}$) @endtex |
---|
1438 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! structural resistance @tex ($s m^{-1}$) @endtex |
---|
1439 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! mean intercellular CO2 concentration |
---|
1440 | !! @tex ($\mu mol mol^{-1}$) @endtex |
---|
1441 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gsmean !! mean stomatal conductance to CO2 (umol m-2 s-1) |
---|
1442 | REAL(r_Std),DIMENSION (kjpindex,nvm), INTENT (out) :: gpp !! Assimilation ((gC m^{-2} dt_sechiba^{-1}), total area) |
---|
1443 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cim !! Intercellular CO2 over nlai |
---|
1444 | ! |
---|
1445 | !! 0.3 Modified variables |
---|
1446 | ! |
---|
1447 | |
---|
1448 | ! |
---|
1449 | !! 0.4 Local variables |
---|
1450 | ! |
---|
1451 | REAL(r_std),DIMENSION (kjpindex) :: vcmax !! maximum rate of carboxylation. Notice that, |
---|
1452 | !! as it is currently written, this is only |
---|
1453 | !! a function of nue and leafN_top |
---|
1454 | !! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex |
---|
1455 | REAL(r_std),DIMENSION (kjpindex,nvm) :: nue !! Nitrogen use Efficiency with impact of leaf age (umol CO2 (gN)-1 s-1) |
---|
1456 | REAL(r_std),DIMENSION (kjpindex,nvm) :: ntot !! Total canopy (eg leaf) N (gN m-2[ground]) |
---|
1457 | REAL(r_std) :: N_profile !! N-content reduction factor for a specific canopy depth (unitless) |
---|
1458 | REAL(r_std),DIMENSION (kjpindex) :: leafN_top !! Leaf N content - top of canopy (gN m-2[leaf]) |
---|
1459 | INTEGER(i_std) :: ji, jv, jl, limit_photo !! indices (unitless) |
---|
1460 | REAL(r_std), DIMENSION(kjpindex,nlai+1) :: info_limitphoto |
---|
1461 | REAL(r_std), DIMENSION(kjpindex,nvm,nlai) :: leaf_ci !! intercellular CO2 concentration (ppm) |
---|
1462 | REAL(r_std), DIMENSION(kjpindex) :: leaf_ci_lowest !! intercellular CO2 concentration at the lowest |
---|
1463 | !! LAI level |
---|
1464 | !! @tex ($\mu mol mol^{-1}$) @endtex |
---|
1465 | INTEGER(i_std), DIMENSION(kjpindex) :: ilai !! counter for loops on LAI levels (unitless) |
---|
1466 | REAL(r_std), DIMENSION(kjpindex) :: zqsvegrap !! relative water quantity in the water |
---|
1467 | !! interception reservoir (0-1,unitless) |
---|
1468 | REAL(r_std) :: speed !! wind speed @tex ($m s^{-1}$) @endtex |
---|
1469 | ! Assimilation |
---|
1470 | LOGICAL, DIMENSION(kjpindex) :: assimilate !! where assimilation is to be calculated |
---|
1471 | !! (unitless) |
---|
1472 | LOGICAL, DIMENSION(kjpindex) :: calculate !! where assimilation is to be calculated for |
---|
1473 | !! in the PFTs loop (unitless) |
---|
1474 | INTEGER(i_std) :: nic,inic,icinic !! counter/indices (unitless) |
---|
1475 | INTEGER(i_std), DIMENSION(kjpindex) :: index_calc !! index (unitless) |
---|
1476 | INTEGER(i_std) :: nia,inia,nina,inina,iainia !! counter/indices (unitless) |
---|
1477 | INTEGER(i_std), DIMENSION(kjpindex) :: index_assi,index_non_assi !! indices (unitless) |
---|
1478 | REAL(r_std), DIMENSION(kjpindex, nlai+1) :: vc2 !! rate of carboxylation (at a specific LAI level) |
---|
1479 | !! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex |
---|
1480 | REAL(r_std), DIMENSION(kjpindex, nlai+1) :: vj2 !! rate of Rubisco regeneration (at a specific LAI |
---|
1481 | !! level) @tex ($\mu mol e- m^{-2} s^{-1}$) @endtex |
---|
1482 | REAL(r_std), DIMENSION(kjpindex, nlai+1) :: assimi !! assimilation (at a specific LAI level) |
---|
1483 | !! @tex ($\mu mol m^{-2} s^{-1}$) @endtex |
---|
1484 | !! (temporary variables) |
---|
1485 | REAL(r_std), DIMENSION(kjpindex) :: gstop !! stomatal conductance to H2O at topmost level |
---|
1486 | !! @tex ($m s^{-1}$) @endtex |
---|
1487 | REAL(r_std), DIMENSION(kjpindex,nlai+1) :: gs !! stomatal conductance to CO2 |
---|
1488 | !! @tex ($\mol m^{-2} s^{-1}$) @endtex |
---|
1489 | REAL(r_std), DIMENSION(kjpindex) :: gamma_star !! CO2 compensation point (ppm) |
---|
1490 | !! @tex ($\mu mol mol^{-1}$) @endtex |
---|
1491 | |
---|
1492 | REAL(r_std), DIMENSION(kjpindex) :: air_relhum !! air relative humidity at 2m |
---|
1493 | !! @tex ($kg kg^{-1}$) @endtex |
---|
1494 | REAL(r_std), DIMENSION(kjpindex) :: VPD !! Vapor Pressure Deficit (kPa) |
---|
1495 | REAL(r_std), DIMENSION(kjpindex) :: water_lim !! water limitation factor (0-1,unitless) |
---|
1496 | |
---|
1497 | REAL(r_std), DIMENSION(kjpindex) :: gstot !! total stomatal conductance to H2O |
---|
1498 | !! Final unit is |
---|
1499 | !! @tex ($m s^{-1}$) @endtex |
---|
1500 | REAL(r_std), DIMENSION(kjpindex) :: assimtot !! total assimilation |
---|
1501 | !! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex |
---|
1502 | REAL(r_std), DIMENSION(kjpindex) :: Rdtot !! Total Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 mâ2 sâ1) |
---|
1503 | REAL(r_std), DIMENSION(kjpindex) :: leaf_gs_top !! leaf stomatal conductance to H2O at topmost level |
---|
1504 | !! @tex ($\mol H2O m^{-2} s^{-1}$) @endtex |
---|
1505 | REAL(r_std), DIMENSION(nlai+1) :: laitab !! tabulated LAI steps @tex ($m^2 m^{-2}$) @endtex |
---|
1506 | REAL(r_std), DIMENSION(kjpindex) :: qsatt !! surface saturated humidity at 2m (??) |
---|
1507 | !! @tex ($g g^{-1}$) @endtex |
---|
1508 | REAL(r_std), DIMENSION(nvm,nlai) :: light !! fraction of light that gets through upper LAI |
---|
1509 | !! levels (0-1,unitless) |
---|
1510 | REAL(r_std), DIMENSION(kjpindex) :: T_Vcmax !! Temperature dependance of Vcmax (unitless) |
---|
1511 | REAL(r_std), DIMENSION(kjpindex) :: S_Vcmax_acclim_temp !! Entropy term for Vcmax |
---|
1512 | !! accounting for acclimation to temperature (J K-1 mol-1) |
---|
1513 | REAL(r_std), DIMENSION(kjpindex) :: T_Jmax !! Temperature dependance of Jmax |
---|
1514 | REAL(r_std), DIMENSION(kjpindex) :: S_Jmax_acclim_temp !! Entropy term for Jmax |
---|
1515 | !! accounting for acclimation toxs temperature (J K-1 mol-1) |
---|
1516 | REAL(r_std), DIMENSION(kjpindex) :: T_gm !! Temperature dependance of gmw |
---|
1517 | REAL(r_std), DIMENSION(kjpindex) :: T_Rd !! Temperature dependance of Rd (unitless) |
---|
1518 | REAL(r_std), DIMENSION(kjpindex) :: T_Kmc !! Temperature dependance of KmC (unitless) |
---|
1519 | REAL(r_std), DIMENSION(kjpindex) :: T_KmO !! Temperature dependance of KmO (unitless) |
---|
1520 | REAL(r_std), DIMENSION(kjpindex) :: T_Sco !! Temperature dependance of Sco |
---|
1521 | REAL(r_std), DIMENSION(kjpindex) :: T_gamma_star !! Temperature dependance of gamma_star (unitless) |
---|
1522 | REAL(r_std), DIMENSION(kjpindex) :: vc !! Maximum rate of Rubisco activity-limited carboxylation (mumol CO2 mâ2 sâ1) |
---|
1523 | REAL(r_std), DIMENSION(kjpindex) :: vj !! Maximum rate of e- transport under saturated light (mumol CO2 mâ2 sâ1) |
---|
1524 | REAL(r_std), DIMENSION(kjpindex) :: gm !! Mesophyll diffusion conductance (molCO2 mâ2 sâ1 barâ1) |
---|
1525 | REAL(r_std), DIMENSION(kjpindex) :: g0var |
---|
1526 | REAL(r_std), DIMENSION(kjpindex,nlai+1) :: Rd !! Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 mâ2 sâ1) |
---|
1527 | REAL(r_std), DIMENSION(kjpindex) :: Kmc !! MichaelisâMenten constant of Rubisco for CO2 (mubar) |
---|
1528 | REAL(r_std), DIMENSION(kjpindex) :: KmO !! MichaelisâMenten constant of Rubisco for O2 (mubar) |
---|
1529 | REAL(r_std), DIMENSION(kjpindex) :: Sco !! Relative CO2 /O2 specificity factor for Rubisco (bar bar-1) |
---|
1530 | REAL(r_std), DIMENSION(kjpindex) :: gb_co2 !! Boundary-layer conductance (molCO2 mâ2 sâ1 barâ1) |
---|
1531 | REAL(r_std), DIMENSION(kjpindex) :: gb_h2o !! Boundary-layer conductance (molH2O mâ2 sâ1 barâ1) |
---|
1532 | REAL(r_std), DIMENSION(kjpindex) :: fvpd !! Factor for describing the effect of leaf-to-air vapour difference on gs (-) |
---|
1533 | REAL(r_std), DIMENSION(kjpindex) :: low_gamma_star !! Half of the reciprocal of Sc/o (bar bar-1) |
---|
1534 | REAL(r_std) :: fcyc !! Fraction of electrons at PSI that follow cyclic transport around PSI (-) |
---|
1535 | REAL(r_std) :: z !! A lumped parameter (see Yin et al. 2009) ( mol mol-1) |
---|
1536 | REAL(r_std) :: Rm !! Day respiration in the mesophyll (umol CO2 mâ2 sâ1) |
---|
1537 | REAL(r_std) :: Cs_star !! Cs -based CO2 compensation point in the absence of Rd (ubar) |
---|
1538 | REAL(r_std), DIMENSION(kjpindex) :: Iabs !! Photon flux density absorbed by leaf photosynthetic pigments (umol photon mâ2 sâ1) |
---|
1539 | REAL(r_std), DIMENSION(kjpindex) :: Jmax !! Maximum value of J under saturated light (umol eâ mâ2 sâ1) |
---|
1540 | REAL(r_std), DIMENSION(kjpindex, nlai+1) :: JJ !! Rate of eâ transport (umol eâ mâ2 sâ1) |
---|
1541 | REAL(r_std) :: J2 !! Rate of all eâ transport through PSII (umol eâ mâ2 sâ1) |
---|
1542 | REAL(r_std) :: VpJ2 !! eâ transport-limited PEP carboxylation rate (umol CO2 mâ2 sâ1) |
---|
1543 | REAL(r_std) :: A_1, A_3 !! Lowest First and third roots of the analytical solution for a general cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 mâ2 sâ1) |
---|
1544 | REAL(r_std) :: A_1_tmp, A_3_tmp !! Temporary First and third roots of the analytical solution for a general cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 mâ2 sâ1) |
---|
1545 | REAL(r_std) :: Obs !! Bundle-sheath oxygen partial pressure (ubar) |
---|
1546 | REAL(r_std), DIMENSION(kjpindex, nlai+1) :: Cc !! Chloroplast CO2 partial pressure (ubar) |
---|
1547 | REAL(r_std) :: ci_star !! Ci -based CO2 compensation point in the absence of Rd (ubar) |
---|
1548 | 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)) |
---|
1549 | REAL(r_std) :: QQ,UU,PSI,x1,x2,x3 !! Variables used for solving the cubic equation (see Yin et al. (2009)) |
---|
1550 | |
---|
1551 | REAL(r_std) :: cresist !! coefficient for resistances (??) |
---|
1552 | REAL(r_std), DIMENSION(kjpindex) :: laisum !! when calculating cim over nlai |
---|
1553 | |
---|
1554 | ! @defgroup Photosynthesis Photosynthesis |
---|
1555 | ! @{ |
---|
1556 | ! 1. Preliminary calculations\n |
---|
1557 | !_ ================================================================================================================================ |
---|
1558 | |
---|
1559 | cim(:,:)=zero |
---|
1560 | leaf_ci(:,:,:) = zero |
---|
1561 | |
---|
1562 | ! |
---|
1563 | ! 1.1 Calculate LAI steps\n |
---|
1564 | ! The integration at the canopy level is done over nlai fixed levels. |
---|
1565 | !! \latexonly |
---|
1566 | !! \input{diffuco_trans_co2_1.1.tex} |
---|
1567 | !! \endlatexonly |
---|
1568 | ! @} |
---|
1569 | ! @codeinc |
---|
1570 | DO jl = 1, nlai+1 |
---|
1571 | laitab(jl) = laimax*(EXP(lai_level_depth*REAL(jl-1,r_std))-1.)/(EXP(lai_level_depth*REAL(nlai,r_std))-un) |
---|
1572 | ENDDO |
---|
1573 | ! @endcodeinc |
---|
1574 | |
---|
1575 | ! @addtogroup Photosynthesis |
---|
1576 | ! @{ |
---|
1577 | ! |
---|
1578 | ! 1.2 Calculate light fraction for each LAI step\n |
---|
1579 | ! The available light follows a simple Beer extinction law. |
---|
1580 | ! The extinction coefficients (ext_coef) are PFT-dependant constants and are defined in constant_co2.f90. |
---|
1581 | !! \latexonly |
---|
1582 | !! \input{diffuco_trans_co2_1.2.tex} |
---|
1583 | !! \endlatexonly |
---|
1584 | ! @} |
---|
1585 | ! @codeinc |
---|
1586 | DO jl = 1, nlai |
---|
1587 | DO jv = 1, nvm |
---|
1588 | light(jv,jl) = exp( -ext_coeff(jv)*laitab(jl) ) |
---|
1589 | ENDDO |
---|
1590 | ENDDO |
---|
1591 | ! @endcodeinc |
---|
1592 | ! |
---|
1593 | ! Photosynthesis parameters |
---|
1594 | ! |
---|
1595 | |
---|
1596 | ! IF (downregulation_co2) THEN |
---|
1597 | ! DO jv= 1, nvm |
---|
1598 | ! vcmax(:,jv) = assim_param(:,jv,ivcmax)*(un-downregulation_co2_coeff(jv)*log(Ca(:)/downregulation_co2_baselevel)) |
---|
1599 | ! ENDDO |
---|
1600 | ! ELSE |
---|
1601 | ! vcmax(:,:) = assim_param(:,:,ivcmax) |
---|
1602 | ! ENDIF |
---|
1603 | |
---|
1604 | |
---|
1605 | nue(:,:) = assim_param(:,:,inue) |
---|
1606 | ntot(:,:) = assim_param(:,:,ileafn) |
---|
1607 | |
---|
1608 | |
---|
1609 | |
---|
1610 | ! DO jv = 1, nvm |
---|
1611 | ! vcmax(:,:) = Vcmax25(jv) |
---|
1612 | ! ENDDO |
---|
1613 | |
---|
1614 | ! @addtogroup Photosynthesis |
---|
1615 | ! @{ |
---|
1616 | ! |
---|
1617 | ! 1.3 Estimate relative humidity of air (for calculation of the stomatal conductance).\n |
---|
1618 | !! \latexonly |
---|
1619 | !! \input{diffuco_trans_co2_1.3.tex} |
---|
1620 | !! \endlatexonly |
---|
1621 | ! @} |
---|
1622 | ! |
---|
1623 | |
---|
1624 | CALL qsatcalc (kjpindex, temp_air, pb, qsatt) |
---|
1625 | air_relhum(:) = & |
---|
1626 | ( qair(:) * pb(:) / (Tetens_1+qair(:)* Tetens_2) ) / & |
---|
1627 | ( qsatt(:)*pb(:) / (Tetens_1+qsatt(:)*Tetens_2 ) ) |
---|
1628 | |
---|
1629 | |
---|
1630 | VPD(:) = ( qsatt(:)*pb(:) / (Tetens_1+qsatt(:)*Tetens_2 ) ) & |
---|
1631 | - ( qair(:) * pb(:) / (Tetens_1+qair(:)* Tetens_2) ) |
---|
1632 | ! VPD is needed in kPa |
---|
1633 | VPD(:) = VPD(:)/10. |
---|
1634 | |
---|
1635 | ! |
---|
1636 | ! 2. beta coefficient for vegetation transpiration |
---|
1637 | ! |
---|
1638 | rstruct(:,1) = rstruct_const(1) |
---|
1639 | rveget(:,:) = undef_sechiba |
---|
1640 | ! |
---|
1641 | vbeta3(:,:) = zero |
---|
1642 | vbeta3pot(:,:) = zero |
---|
1643 | gsmean(:,:) = zero |
---|
1644 | gpp(:,:) = zero |
---|
1645 | ! |
---|
1646 | cimean(:,1) = Ca(:) |
---|
1647 | ! |
---|
1648 | ! @addtogroup Photosynthesis |
---|
1649 | ! @{ |
---|
1650 | ! 2. Loop over vegetation types\n |
---|
1651 | ! @} |
---|
1652 | ! |
---|
1653 | DO jv = 2,nvm |
---|
1654 | gamma_star(:)=zero |
---|
1655 | Kmo(:)=zero |
---|
1656 | Kmc(:)=zero |
---|
1657 | gm(:)=zero |
---|
1658 | g0var(:) =zero |
---|
1659 | |
---|
1660 | Cc(:,:)=zero |
---|
1661 | Vc2(:,:)=zero |
---|
1662 | JJ(:,:)=zero |
---|
1663 | info_limitphoto(:,:)=zero |
---|
1664 | gs(:,:)=zero |
---|
1665 | assimi(:,:)=zero |
---|
1666 | Rd(:,:)=zero |
---|
1667 | |
---|
1668 | ! |
---|
1669 | ! @addtogroup Photosynthesis |
---|
1670 | ! @{ |
---|
1671 | ! |
---|
1672 | ! 2.1 Initializations\n |
---|
1673 | !! \latexonly |
---|
1674 | !! \input{diffuco_trans_co2_2.1.tex} |
---|
1675 | !! \endlatexonly |
---|
1676 | ! @} |
---|
1677 | ! |
---|
1678 | ! beta coefficient for vegetation transpiration |
---|
1679 | ! |
---|
1680 | rstruct(:,jv) = rstruct_const(jv) |
---|
1681 | cimean(:,jv) = Ca(:) |
---|
1682 | |
---|
1683 | WHERE (lai(:,jv) .GT. min_sechiba) |
---|
1684 | leafN_top(:) = ntot(:,jv) * ext_coeff_N(jv) / ( 1. - exp(-ext_coeff_N(jv) * lai(:,jv)) ) |
---|
1685 | ELSEWHERE |
---|
1686 | leafN_top(:) = zero |
---|
1687 | ENDWHERE |
---|
1688 | |
---|
1689 | vcmax(:) = nue(:,jv)*leafN_top(:) |
---|
1690 | |
---|
1691 | ! |
---|
1692 | !! mask that contains points where there is photosynthesis |
---|
1693 | !! For the sake of vectorisation [DISPENSABLE], computations are done only for convenient points. |
---|
1694 | !! nia is the number of points where the assimilation is calculated and nina the number of points where photosynthesis is not |
---|
1695 | !! calculated (based on criteria on minimum or maximum values on LAI, vegetation fraction, shortwave incoming radiation, |
---|
1696 | !! temperature and relative humidity). |
---|
1697 | !! For the points where assimilation is not calculated, variables are initialized to specific values. |
---|
1698 | !! The assimilate(kjpindex) array contains the logical value (TRUE/FALSE) relative to this photosynthesis calculation. |
---|
1699 | !! The index_assi(kjpindex) array indexes the nia points with assimilation, whereas the index_no_assi(kjpindex) array indexes |
---|
1700 | !! the nina points with no assimilation. |
---|
1701 | nia=0 |
---|
1702 | nina=0 |
---|
1703 | ! |
---|
1704 | DO ji=1,kjpindex |
---|
1705 | ! |
---|
1706 | IF ( ( lai(ji,jv) .GT. 0.01 ) .AND. & |
---|
1707 | ( veget_max(ji,jv) .GT. min_sechiba ) ) THEN |
---|
1708 | |
---|
1709 | IF ( ( veget(ji,jv) .GT. min_sechiba ) .AND. & |
---|
1710 | ( swdown(ji) .GT. min_sechiba ) .AND. & |
---|
1711 | ( humrel(ji,jv) .GT. min_sechiba) .AND. & |
---|
1712 | ( temp_growth(ji) .GT. tphoto_min(jv) ) .AND. & |
---|
1713 | ( temp_growth(ji) .LT. tphoto_max(jv) ) ) THEN |
---|
1714 | ! |
---|
1715 | assimilate(ji) = .TRUE. |
---|
1716 | nia=nia+1 |
---|
1717 | index_assi(nia)=ji |
---|
1718 | ! |
---|
1719 | ELSE |
---|
1720 | ! |
---|
1721 | assimilate(ji) = .FALSE. |
---|
1722 | nina=nina+1 |
---|
1723 | index_non_assi(nina)=ji |
---|
1724 | ! |
---|
1725 | ENDIF |
---|
1726 | ELSE |
---|
1727 | ! |
---|
1728 | assimilate(ji) = .FALSE. |
---|
1729 | nina=nina+1 |
---|
1730 | index_non_assi(nina)=ji |
---|
1731 | ! |
---|
1732 | ENDIF |
---|
1733 | ! |
---|
1734 | |
---|
1735 | ENDDO |
---|
1736 | ! |
---|
1737 | |
---|
1738 | gstot(:) = zero |
---|
1739 | gstop(:) = zero |
---|
1740 | assimtot(:) = zero |
---|
1741 | Rdtot(:)=zero |
---|
1742 | leaf_gs_top(:) = zero |
---|
1743 | ! |
---|
1744 | zqsvegrap(:) = zero |
---|
1745 | WHERE (qsintmax(:,jv) .GT. min_sechiba) |
---|
1746 | !! relative water quantity in the water interception reservoir |
---|
1747 | zqsvegrap(:) = MAX(zero, qsintveg(:,jv) / qsintmax(:,jv)) |
---|
1748 | ENDWHERE |
---|
1749 | ! |
---|
1750 | !! Calculates the water limitation factor. |
---|
1751 | water_lim(:) = humrel(:,jv) |
---|
1752 | |
---|
1753 | ! give a default value of ci for all pixel that do not assimilate |
---|
1754 | DO jl=1,nlai |
---|
1755 | DO inina=1,nina |
---|
1756 | leaf_ci(index_non_assi(inina),jv,jl) = Ca(index_non_assi(inina)) |
---|
1757 | ENDDO |
---|
1758 | ENDDO |
---|
1759 | ! |
---|
1760 | ilai(:) = 1 |
---|
1761 | ! |
---|
1762 | ! Here is the calculation of assimilation and stomatal conductance |
---|
1763 | ! based on the work of Farquahr, von Caemmerer and Berry (FvCB model) |
---|
1764 | ! as described in Yin et al. 2009 |
---|
1765 | ! Yin et al. developed a extended version of the FvCB model for C4 plants |
---|
1766 | ! and proposed an analytical solution for both photosynthesis pathways (C3 and C4) |
---|
1767 | ! Photosynthetic parameters used are those reported in Yin et al. |
---|
1768 | ! Except For Vcmax25, relationships between Vcmax25 and Jmax25 for which we use |
---|
1769 | ! Medlyn et al. (2002) and Kattge & Knorr (2007) |
---|
1770 | ! Because these 2 references do not consider mesophyll conductance, we neglect this term |
---|
1771 | ! in the formulations developed by Yin et al. |
---|
1772 | ! Consequently, gm (the mesophyll conductance) tends to the infinite |
---|
1773 | ! This is of importance because as stated by Kattge & Knorr and Medlyn et al., |
---|
1774 | ! values of Vcmax and Jmax derived with different model parametrizations are not |
---|
1775 | ! directly comparable and the published values of Vcmax and Jmax had to be standardized |
---|
1776 | ! to one consistent formulation and parametrization |
---|
1777 | |
---|
1778 | ! See eq. 6 of Yin et al. (2009) |
---|
1779 | ! Parametrization of Medlyn et al. (2002) - from Bernacchi et al. (2001) |
---|
1780 | T_KmC(:) = Arrhenius(kjpindex,temp_air,298.,E_KmC(jv)) |
---|
1781 | T_KmO(:) = Arrhenius(kjpindex,temp_air,298.,E_KmO(jv)) |
---|
1782 | T_Sco(:) = Arrhenius(kjpindex,temp_air,298.,E_Sco(jv)) |
---|
1783 | T_gamma_star(:) = Arrhenius(kjpindex,temp_air,298.,E_gamma_star(jv)) |
---|
1784 | |
---|
1785 | |
---|
1786 | ! Parametrization of Yin et al. (2009) - from Bernacchi et al. (2001) |
---|
1787 | T_Rd(:) = Arrhenius(kjpindex,temp_air,298.,E_Rd(jv)) |
---|
1788 | |
---|
1789 | |
---|
1790 | ! For C3 plants, we assume that the Entropy term for Vcmax and Jmax |
---|
1791 | ! acclimates to temperature as shown by Kattge & Knorr (2007) - Eq. 9 and 10 |
---|
1792 | ! and that Jmax and Vcmax respond to temperature following a modified Arrhenius function |
---|
1793 | ! (with a decrease of these parameters for high temperature) as in Medlyn et al. (2002) |
---|
1794 | ! and Kattge & Knorr (2007). |
---|
1795 | ! In Yin et al. (2009), temperature dependance to Vcmax is based only on a Arrhenius function |
---|
1796 | ! Concerning this apparent unconsistency, have a look to the section 'Limitation of |
---|
1797 | ! Photosynthesis by gm' of Bernacchi (2002) that may provide an explanation |
---|
1798 | |
---|
1799 | ! Growth temperature tested by Kattge & Knorr range from 11 to 35°C |
---|
1800 | ! So, we limit the relationship between these lower and upper limits |
---|
1801 | S_Jmax_acclim_temp(:) = aSJ(jv) + bSJ(jv) * MAX(11., MIN(temp_growth(:),35.)) |
---|
1802 | T_Jmax(:) = Arrhenius_modified(kjpindex,temp_air,298.,E_Jmax(jv),D_Jmax(jv),S_Jmax_acclim_temp) |
---|
1803 | |
---|
1804 | S_Vcmax_acclim_temp(:) = aSV(jv) + bSV(jv) * MAX(11., MIN(temp_growth(:),35.)) |
---|
1805 | T_Vcmax(:) = Arrhenius_modified(kjpindex,temp_air,298.,E_Vcmax(jv),D_Vcmax(jv),S_Vcmax_acclim_temp) |
---|
1806 | |
---|
1807 | |
---|
1808 | |
---|
1809 | vc(:) = vcmax(:) * T_Vcmax(:) |
---|
1810 | ! As shown by Kattge & Knorr (2007), we make use |
---|
1811 | ! of Jmax25/Vcmax25 ratio (rJV) that acclimates to temperature for C3 plants |
---|
1812 | ! rJV is written as a function of the growth temperature |
---|
1813 | ! rJV = arJV + brJV * T_month |
---|
1814 | ! See eq. 10 of Kattge & Knorr (2007) |
---|
1815 | ! and Table 3 for Values of arJV anf brJV |
---|
1816 | ! Growth temperature is monthly temperature (expressed in °C) - See first paragraph of |
---|
1817 | ! section Methods/Data of Kattge & Knorr |
---|
1818 | vj(:) = ( arJV(jv) + brJV(jv) * MAX(11., MIN(temp_growth(:),35.)) ) * vcmax(:) * T_Jmax(:) |
---|
1819 | |
---|
1820 | T_gm(:) = Arrhenius_modified(kjpindex,temp_air,298.,E_gm(jv),D_gm(jv),S_gm(jv)) |
---|
1821 | gm(:) = gm25(jv) * T_gm(:) * MAX(1-stress_gm(jv), water_lim(:)) |
---|
1822 | |
---|
1823 | g0var(:) = g0(jv)* MAX(1-stress_gs(jv), water_lim(:)) |
---|
1824 | ! @endcodeinc |
---|
1825 | ! |
---|
1826 | KmC(:)=KmC25(jv)*T_KmC(:) |
---|
1827 | KmO(:)=KmO25(jv)*T_KmO(:) |
---|
1828 | Sco(:)=Sco25(jv)*T_sco(:) |
---|
1829 | gamma_star(:) = gamma_star25(jv)*T_gamma_star(:) |
---|
1830 | |
---|
1831 | |
---|
1832 | |
---|
1833 | ! low_gamma_star is defined by Yin et al. (2009) |
---|
1834 | ! as the half of the reciprocal of Sco - See Table 2 |
---|
1835 | low_gamma_star(:) = 0.5 / Sco(:) |
---|
1836 | |
---|
1837 | ! VPD expressed in kPa |
---|
1838 | ! Note : MIN(1.-min_sechiba,MAX(min_sechiba,(a1(jv) - b1(jv) * VPD(:)))) is always between 0-1 not including 0 and 1 |
---|
1839 | fvpd(:) = 1. / ( 1. / MIN(1.-min_sechiba,MAX(min_sechiba,(a1(jv) - b1(jv) * VPD(:)))) - 1. ) & |
---|
1840 | * MAX(1-stress_gs(jv), water_lim(:)) |
---|
1841 | |
---|
1842 | ! leaf boundary layer conductance |
---|
1843 | ! conversion from a conductance in (m s-1) to (mol H2O m-2 s-1) |
---|
1844 | ! from Pearcy et al. (1991, see below) |
---|
1845 | gb_h2o(:) = gb_ref * 44.6 * (tp_00/temp_air(:)) * (pb(:)/pb_std) |
---|
1846 | |
---|
1847 | ! conversion from (mol H2O m-2 s-1) to (mol CO2 m-2 s-1) |
---|
1848 | gb_co2(:) = gb_h2o(:) / ratio_H2O_to_CO2 |
---|
1849 | |
---|
1850 | ! |
---|
1851 | ! @addtogroup Photosynthesis |
---|
1852 | ! @{ |
---|
1853 | ! |
---|
1854 | ! 2.4 Loop over LAI discretized levels to estimate assimilation and conductance\n |
---|
1855 | ! @} |
---|
1856 | ! |
---|
1857 | !! The calculate(kjpindex) array is of type logical to indicate wether we have to sum over this LAI fixed level (the LAI of |
---|
1858 | !! the point for the PFT is lower or equal to the LAI level value). The number of such points is incremented in nic and the |
---|
1859 | !! corresponding point is indexed in the index_calc array. |
---|
1860 | JJ(:,:)=zero |
---|
1861 | vc2(:,:)=zero |
---|
1862 | vj2(:,:)=zero |
---|
1863 | Cc(:,:)=zero |
---|
1864 | gs(:,:)=zero |
---|
1865 | assimi(:,:)=zero |
---|
1866 | Rd(:,:)=zero |
---|
1867 | |
---|
1868 | DO jl = 1, nlai |
---|
1869 | ! |
---|
1870 | nic=0 |
---|
1871 | calculate(:) = .FALSE. |
---|
1872 | ! |
---|
1873 | IF (nia .GT. 0) then |
---|
1874 | DO inia=1,nia |
---|
1875 | calculate(index_assi(inia)) = (laitab(jl) .LE. lai(index_assi(inia),jv) ) |
---|
1876 | IF ( calculate(index_assi(inia)) ) THEN |
---|
1877 | nic=nic+1 |
---|
1878 | index_calc(nic)=index_assi(inia) |
---|
1879 | ENDIF |
---|
1880 | ENDDO |
---|
1881 | ENDIF |
---|
1882 | ! |
---|
1883 | ! @addtogroup Photosynthesis |
---|
1884 | ! @{ |
---|
1885 | ! |
---|
1886 | ! 2.4.1 Vmax is scaled into the canopy due to reduction of nitrogen |
---|
1887 | !! (Johnson and Thornley,1984).\n |
---|
1888 | !! \latexonly |
---|
1889 | !! \input{diffuco_trans_co2_2.4.1.tex} |
---|
1890 | !! \endlatexonly |
---|
1891 | ! @} |
---|
1892 | ! |
---|
1893 | N_profile = exp( -ext_coeff_N(jv)*laitab(jl) ) |
---|
1894 | ! |
---|
1895 | |
---|
1896 | vc2(:,jl) = vc(:) * N_profile * MAX(1-stress_vcmax(jv),water_lim(:)) |
---|
1897 | vj2(:,jl) = vj(:) * N_profile * MAX(1-stress_vcmax(jv),water_lim(:)) |
---|
1898 | |
---|
1899 | ! see Comment in legend of Fig. 6 of Yin et al. (2009) |
---|
1900 | ! Rd25 is assumed to equal 0.01 Vcmax25 |
---|
1901 | Rd(:,jl) = vcmax(:) * N_profile * 0.01 * T_Rd(:) * MAX(1-stress_vcmax(jv), water_lim(:)) |
---|
1902 | |
---|
1903 | Iabs(:)=swdown(:)*W_to_mol*RG_to_PAR*ext_coeff(jv)*light(jv,jl) |
---|
1904 | |
---|
1905 | ! eq. 4 of Yin et al (2009) |
---|
1906 | Jmax(:)=vj2(:,jl) |
---|
1907 | JJ(:,jl) = ( alpha_LL(jv) * Iabs(:) + Jmax(:) - sqrt((alpha_LL(jv) * Iabs(:) + Jmax(:) )**2. & |
---|
1908 | - 4 * theta(jv) * Jmax(:) * alpha_LL(jv) * Iabs(:)) ) & |
---|
1909 | / ( 2 * theta(jv)) |
---|
1910 | |
---|
1911 | ! |
---|
1912 | IF ( is_c4(jv) ) THEN |
---|
1913 | ! |
---|
1914 | ! @addtogroup Photosynthesis |
---|
1915 | ! @{ |
---|
1916 | ! |
---|
1917 | ! 2.4.2 Assimilation for C4 plants (Collatz et al., 1992)\n |
---|
1918 | !! \latexonly |
---|
1919 | !! \input{diffuco_trans_co2_2.4.2.tex} |
---|
1920 | !! \endlatexonly |
---|
1921 | ! @} |
---|
1922 | ! |
---|
1923 | ! |
---|
1924 | ! |
---|
1925 | IF (nic .GT. 0) THEN |
---|
1926 | DO inic=1,nic |
---|
1927 | |
---|
1928 | ! Analytical resolution of the Assimilation based Yin et al. (2009) |
---|
1929 | icinic=index_calc(inic) |
---|
1930 | |
---|
1931 | ! Eq. 28 of Yin et al. (2009) |
---|
1932 | fcyc= 1. - ( 4.*(1.-fpsir(jv))*(1.+fQ(jv)) + 3.*h_protons(jv)*fpseudo(jv) ) / & |
---|
1933 | ( 3.*h_protons(jv) - 4.*(1.-fpsir(jv))) |
---|
1934 | |
---|
1935 | ! See paragraph after eq. (20b) of Yin et al. |
---|
1936 | Rm=Rd(icinic,jl)/2. |
---|
1937 | |
---|
1938 | ! We assume that cs_star equals ci_star (see Comment in legend of Fig. 6 of Yin et al. (2009) |
---|
1939 | ! Equation 26 of Yin et al. (2009) |
---|
1940 | Cs_star = (gbs(jv) * low_gamma_star(icinic) * Oi - & |
---|
1941 | ( 1. + low_gamma_star(icinic) * alpha(jv) / 0.047) * Rd(icinic,jl) + Rm ) & |
---|
1942 | / ( gbs(jv) + kp(jv) ) |
---|
1943 | |
---|
1944 | ! eq. 11 of Yin et al (2009) |
---|
1945 | J2 = JJ(icinic,jl) / ( 1. - fpseudo(jv) / ( 1. - fcyc ) ) |
---|
1946 | |
---|
1947 | ! Equation right after eq. (20d) of Yin et al. (2009) |
---|
1948 | z = ( 2. + fQ(jv) - fcyc ) / ( h_protons(jv) * (1. - fcyc )) |
---|
1949 | |
---|
1950 | VpJ2 = fpsir(jv) * J2 * z / 2. |
---|
1951 | |
---|
1952 | A_3=9999. |
---|
1953 | |
---|
1954 | ! See eq. right after eq. 18 of Yin et al. (2009) |
---|
1955 | DO limit_photo=1,2 |
---|
1956 | ! Is Vc limiting the Assimilation |
---|
1957 | IF ( limit_photo .EQ. 1 ) THEN |
---|
1958 | a = 1. + kp(jv) / gbs(jv) |
---|
1959 | b = 0. |
---|
1960 | x1 = Vc2(icinic,jl) |
---|
1961 | x2 = KmC(icinic)/KmO(icinic) |
---|
1962 | x3 = KmC(icinic) |
---|
1963 | ! Is J limiting the Assimilation |
---|
1964 | ELSE |
---|
1965 | a = 1. |
---|
1966 | b = VpJ2 |
---|
1967 | x1 = (1.- fpsir(jv)) * J2 * z / 3. |
---|
1968 | x2 = 7. * low_gamma_star(icinic) / 3. |
---|
1969 | x3 = 0. |
---|
1970 | ENDIF |
---|
1971 | |
---|
1972 | m=fvpd(icinic)-g0var(icinic)/gb_co2(icinic) |
---|
1973 | d=g0var(icinic)*(Ca(icinic)-Cs_star) + fvpd(icinic)*Rd(icinic,jl) |
---|
1974 | f=(b-Rm-low_gamma_star(icinic)*Oi*gbs(jv))*x1*d + a*gbs(jv)*x1*Ca(icinic)*d |
---|
1975 | j=(b-Rm+gbs(jv)*x3 + x2*gbs(jv)*Oi)*m + (alpha(jv)*x2/0.047-1.)*d & |
---|
1976 | + a*gbs(jv)*(Ca(icinic)*m - d/gb_co2(icinic) - (Ca(icinic) - Cs_star )) |
---|
1977 | |
---|
1978 | g=(b-Rm-low_gamma_star(icinic)*Oi*gbs(jv))*x1*m - (alpha(jv)*low_gamma_star(icinic)/0.047+1.)*x1*d & |
---|
1979 | + a*gbs(jv)*x1*(Ca(icinic)*m - d/gb_co2(icinic) - (Ca(icinic)-Cs_star )) |
---|
1980 | |
---|
1981 | h=-((alpha(jv)*low_gamma_star(icinic)/0.047+1.)*x1*m + (a*gbs(jv)*x1*(m-1.))/gb_co2(icinic) ) |
---|
1982 | i= ( b-Rm + gbs(jv)*x3 + x2*gbs(jv)*Oi )*d + a*gbs(jv)*Ca(icinic)*d |
---|
1983 | l= ( alpha(jv)*x2/0.047 - 1.)*m - (a*gbs(jv)*(m-1.))/gb_co2(icinic) |
---|
1984 | |
---|
1985 | p = (j-(h-l*Rd(icinic,jl))) / l |
---|
1986 | q = (i+j*Rd(icinic,jl)-g) / l |
---|
1987 | r = -(f-i*Rd(icinic,jl)) / l |
---|
1988 | |
---|
1989 | ! See Yin et al. (2009) and Baldocchi (1994) |
---|
1990 | QQ = ( (p**2._r_std) - 3._r_std * q) / 9._r_std |
---|
1991 | UU = ( 2._r_std* (p**3._r_std) - 9._r_std *p*q + 27._r_std *r) /54._r_std |
---|
1992 | |
---|
1993 | IF ( (QQ .GE. 0._r_std) .AND. (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) ) THEN |
---|
1994 | PSI = ACOS(UU/(QQ**1.5_r_std)) |
---|
1995 | A_3_tmp = -2._r_std * SQRT(QQ) * COS(( PSI + 4._r_std * PI)/3._r_std ) - p / 3._r_std |
---|
1996 | IF (( A_3_tmp .LT. A_3 )) THEN |
---|
1997 | A_3 = A_3_tmp |
---|
1998 | info_limitphoto(icinic,jl)=2. |
---|
1999 | ELSE |
---|
2000 | ! In case, J is not limiting the assimilation |
---|
2001 | ! we have to re-initialise a, b, x1, x2 and x3 values |
---|
2002 | ! in agreement with a Vc-limited assimilation |
---|
2003 | a = 1. + kp(jv) / gbs(jv) |
---|
2004 | b = 0. |
---|
2005 | x1 = Vc2(icinic,jl) |
---|
2006 | x2 = KmC(icinic)/KmO(icinic) |
---|
2007 | x3 = KmC(icinic) |
---|
2008 | info_limitphoto(icinic,jl)=1. |
---|
2009 | ENDIF |
---|
2010 | ENDIF |
---|
2011 | |
---|
2012 | IF ( ( A_3 .EQ. 9999. ) .OR. ( A_3 .LT. (-Rd(icinic,jl)) ) ) THEN |
---|
2013 | IF ( printlev>=4 ) THEN |
---|
2014 | WRITE(numout,*) 'We have a problem in diffuco_trans_co2 for A_3' |
---|
2015 | WRITE(numout,*) 'no real positive solution found for pft:',jv |
---|
2016 | WRITE(numout,*) 'temp_air:',temp_air(icinic) |
---|
2017 | WRITE(numout,*) 'vpd:',vpd(icinic) |
---|
2018 | END IF |
---|
2019 | A_3 = -Rd(icinic,jl) |
---|
2020 | ENDIF |
---|
2021 | assimi(icinic,jl) = A_3 |
---|
2022 | |
---|
2023 | IF ( ABS( assimi(icinic,jl) + Rd(icinic,jl) ) .LT. min_sechiba ) THEN |
---|
2024 | gs(icinic,jl) = g0var(icinic) |
---|
2025 | !leaf_ci keeps its initial value (Ca). |
---|
2026 | ELSE |
---|
2027 | ! Eq. 24 of Yin et al. (2009) |
---|
2028 | Obs = ( alpha(jv) * assimi(icinic,jl) ) / ( 0.047 * gbs(jv) ) + Oi |
---|
2029 | ! Eq. 23 of Yin et al. (2009) |
---|
2030 | Cc(icinic,jl) = ( ( assimi(icinic,jl) + Rd(icinic,jl) ) * ( x2 * Obs + x3 ) + low_gamma_star(icinic) & |
---|
2031 | * Obs * x1 ) & |
---|
2032 | / MAX(min_sechiba, x1 - ( assimi(icinic,jl) + Rd(icinic,jl) )) |
---|
2033 | ! Eq. 22 of Yin et al. (2009) |
---|
2034 | leaf_ci(icinic,jv,jl) = ( Cc(icinic,jl) - ( b - assimi(icinic,jl) - Rm ) / gbs(jv) ) / a |
---|
2035 | ! Eq. 25 of Yin et al. (2009) |
---|
2036 | ! It should be Cs instead of Ca but it seems that |
---|
2037 | ! other equations in Appendix C make use of Ca |
---|
2038 | gs(icinic,jl) = g0var(icinic) + ( assimi(icinic,jl) + Rd(icinic,jl) ) / & |
---|
2039 | ( Ca(icinic) - Cs_star ) * fvpd(icinic) |
---|
2040 | ENDIF |
---|
2041 | ENDDO |
---|
2042 | ENDDO |
---|
2043 | ENDIF |
---|
2044 | ELSE |
---|
2045 | ! |
---|
2046 | ! @addtogroup Photosynthesis |
---|
2047 | ! @{ |
---|
2048 | ! |
---|
2049 | ! 2.4.3 Assimilation for C3 plants (Farqhuar et al., 1980)\n |
---|
2050 | !! \latexonly |
---|
2051 | !! \input{diffuco_trans_co2_2.4.3.tex} |
---|
2052 | !! \endlatexonly |
---|
2053 | ! @} |
---|
2054 | ! |
---|
2055 | ! |
---|
2056 | IF (nic .GT. 0) THEN |
---|
2057 | DO inic=1,nic |
---|
2058 | icinic=index_calc(inic) |
---|
2059 | |
---|
2060 | A_1=9999. |
---|
2061 | |
---|
2062 | ! See eq. right after eq. 18 of Yin et al. (2009) |
---|
2063 | DO limit_photo=1,2 |
---|
2064 | ! Is Vc limiting the Assimilation |
---|
2065 | IF ( limit_photo .EQ. 1 ) THEN |
---|
2066 | x1 = vc2(icinic,jl) |
---|
2067 | ! It should be O not Oi (comment from Vuichard) |
---|
2068 | x2 = KmC(icinic) * ( 1. + 2*gamma_star(icinic)*Sco(icinic) / KmO(icinic) ) |
---|
2069 | ! Is J limiting the Assimilation |
---|
2070 | ELSE |
---|
2071 | x1 = JJ(icinic,jl)/4. |
---|
2072 | x2 = 2. * gamma_star(icinic) |
---|
2073 | ENDIF |
---|
2074 | |
---|
2075 | |
---|
2076 | ! See Appendix B of Yin et al. (2009) |
---|
2077 | a = g0var(icinic) * ( x2 + gamma_star(icinic) ) + & |
---|
2078 | ( g0var(icinic) / gm(icinic) + fvpd(icinic) ) * ( x1 - Rd(icinic,jl) ) |
---|
2079 | b = Ca(icinic) * ( x1 - Rd(icinic,jl) ) - gamma_star(icinic) * x1 - Rd(icinic,jl) * x2 |
---|
2080 | c = Ca(icinic) + x2 + ( 1./gm(icinic) + 1./gb_co2(icinic) ) * ( x1 - Rd(icinic,jl) ) |
---|
2081 | d = x2 + gamma_star(icinic) + ( x1 - Rd(icinic,jl) ) / gm(icinic) |
---|
2082 | m = 1./gm(icinic) + ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * ( 1./gm(icinic) + 1./gb_co2(icinic) ) |
---|
2083 | |
---|
2084 | p = - ( d + (x1 - Rd(icinic,jl) ) / gm(icinic) + a * ( 1./gm(icinic) + 1./gb_co2(icinic) ) + & |
---|
2085 | ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * c ) / m |
---|
2086 | |
---|
2087 | q = ( d * ( x1 - Rd(icinic,jl) ) + a*c + ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * b ) / m |
---|
2088 | r = - a * b / m |
---|
2089 | |
---|
2090 | ! See Yin et al. (2009) |
---|
2091 | QQ = ( (p**2._r_std) - 3._r_std * q) / 9._r_std |
---|
2092 | UU = ( 2._r_std* (p**3._r_std) - 9._r_std *p*q + 27._r_std *r) /54._r_std |
---|
2093 | |
---|
2094 | IF ( (QQ .GE. 0._r_std) .AND. (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) ) THEN |
---|
2095 | PSI = ACOS(UU/(QQ**1.5_r_std)) |
---|
2096 | A_1_tmp = -2._r_std * SQRT(QQ) * COS( PSI / 3._r_std ) - p / 3._r_std |
---|
2097 | IF (( A_1_tmp .LT. A_1 )) THEN |
---|
2098 | A_1 = A_1_tmp |
---|
2099 | info_limitphoto(icinic,jl)=2. |
---|
2100 | ELSE |
---|
2101 | ! In case, J is not limiting the assimilation |
---|
2102 | ! we have to re-initialise x1 and x2 values |
---|
2103 | ! in agreement with a Vc-limited assimilation |
---|
2104 | x1 = vc2(icinic,jl) |
---|
2105 | ! It should be O not Oi (comment from Vuichard) |
---|
2106 | x2 = KmC(icinic) * ( 1. + 2*gamma_star(icinic)*Sco(icinic) / KmO(icinic) ) |
---|
2107 | info_limitphoto(icinic,jl)=1. |
---|
2108 | ENDIF |
---|
2109 | ENDIF |
---|
2110 | ENDDO |
---|
2111 | IF ( (A_1 .EQ. 9999.) .OR. ( A_1 .LT. (-Rd(icinic,jl)) ) ) THEN |
---|
2112 | IF ( printlev>=4 ) THEN |
---|
2113 | WRITE(numout,*) 'We have a problem in diffuco_trans_co2 for A_1' |
---|
2114 | WRITE(numout,*) 'no real positive solution found for pft:',jv |
---|
2115 | WRITE(numout,*) 'temp_air:',temp_air(icinic) |
---|
2116 | WRITE(numout,*) 'vpd:',vpd(icinic) |
---|
2117 | END IF |
---|
2118 | A_1 = -Rd(icinic,jl) |
---|
2119 | ENDIF |
---|
2120 | assimi(icinic,jl) = A_1 |
---|
2121 | |
---|
2122 | IF ( ABS( assimi(icinic,jl) + Rd(icinic,jl) ) .LT. min_sechiba ) THEN |
---|
2123 | gs(icinic,jl) = g0var(icinic) |
---|
2124 | ELSE |
---|
2125 | ! Eq. 18 of Yin et al. (2009) |
---|
2126 | Cc(icinic,jl) = ( gamma_star(icinic) * x1 + ( assimi(icinic,jl) + Rd(icinic,jl) ) * x2 ) & |
---|
2127 | / MAX( min_sechiba, x1 - ( assimi(icinic,jl) + Rd(icinic,jl) ) ) |
---|
2128 | ! Eq. 17 of Yin et al. (2009) |
---|
2129 | leaf_ci(icinic,jv,jl) = Cc(icinic,jl) + assimi(icinic,jl) / gm(icinic) |
---|
2130 | ! See eq. right after eq. 15 of Yin et al. (2009) |
---|
2131 | ci_star = gamma_star(icinic) - Rd(icinic,jl) / gm(icinic) |
---|
2132 | ! |
---|
2133 | ! Eq. 15 of Yin et al. (2009) |
---|
2134 | gs(icinic,jl) = g0var(icinic) + ( assimi(icinic,jl) + Rd(icinic,jl) ) / ( leaf_ci(icinic,jv,jl) & |
---|
2135 | - ci_star ) * fvpd(icinic) |
---|
2136 | ENDIF |
---|
2137 | ENDDO |
---|
2138 | ENDIF |
---|
2139 | ENDIF |
---|
2140 | ! |
---|
2141 | IF (nic .GT. 0) THEN |
---|
2142 | ! |
---|
2143 | DO inic=1,nic |
---|
2144 | ! |
---|
2145 | ! @addtogroup Photosynthesis |
---|
2146 | ! @{ |
---|
2147 | ! |
---|
2148 | !! 2.4.4 Estimatation of the stomatal conductance (Ball et al., 1987).\n |
---|
2149 | !! \latexonly |
---|
2150 | !! \input{diffuco_trans_co2_2.4.4.tex} |
---|
2151 | !! \endlatexonly |
---|
2152 | ! @} |
---|
2153 | ! |
---|
2154 | icinic=index_calc(inic) |
---|
2155 | ! |
---|
2156 | ! keep stomatal conductance of topmost level |
---|
2157 | ! |
---|
2158 | IF ( jl .EQ. 1 ) THEN |
---|
2159 | leaf_gs_top(icinic) = gs(icinic,jl) |
---|
2160 | ! |
---|
2161 | ENDIF |
---|
2162 | ! |
---|
2163 | ! @addtogroup Photosynthesis |
---|
2164 | ! @{ |
---|
2165 | ! |
---|
2166 | !! 2.4.5 Integration at the canopy level\n |
---|
2167 | !! \latexonly |
---|
2168 | !! \input{diffuco_trans_co2_2.4.5.tex} |
---|
2169 | !! \endlatexonly |
---|
2170 | ! @} |
---|
2171 | ! total assimilation and conductance |
---|
2172 | assimtot(icinic) = assimtot(icinic) + & |
---|
2173 | assimi(icinic,jl) * (laitab(jl+1)-laitab(jl)) |
---|
2174 | Rdtot(icinic) = Rdtot(icinic) + & |
---|
2175 | Rd(icinic,jl) * (laitab(jl+1)-laitab(jl)) |
---|
2176 | gstot(icinic) = gstot(icinic) + & |
---|
2177 | gs(icinic,jl) * (laitab(jl+1)-laitab(jl)) |
---|
2178 | ! |
---|
2179 | ilai(icinic) = jl |
---|
2180 | ! |
---|
2181 | ENDDO |
---|
2182 | ! |
---|
2183 | ENDIF |
---|
2184 | ENDDO ! loop over LAI steps |
---|
2185 | |
---|
2186 | IF(jv==testpft) THEN |
---|
2187 | CALL histwrite_p(hist_id, 'Cc', kjit, Cc, kjpindex*(nlai+1), indexlai) |
---|
2188 | CALL histwrite_p(hist_id, 'Vc', kjit, Vc2, kjpindex*(nlai+1), indexlai) |
---|
2189 | CALL histwrite_p(hist_id, 'Vj', kjit, JJ, kjpindex*(nlai+1), indexlai) |
---|
2190 | CALL histwrite_p(hist_id, 'limitphoto', kjit, info_limitphoto, kjpindex*(nlai+1), indexlai) |
---|
2191 | CALL histwrite_p(hist_id, 'gammastar', kjit, gamma_star, kjpindex,index) |
---|
2192 | CALL histwrite_p(hist_id, 'Kmo', kjit, Kmo, kjpindex,index) |
---|
2193 | CALL histwrite_p(hist_id, 'Kmc', kjit, Kmc, kjpindex,index) |
---|
2194 | CALL histwrite_p(hist_id, 'gm', kjit, gm, kjpindex, index) |
---|
2195 | CALL histwrite_p(hist_id, 'gs', kjit, gs, kjpindex*(nlai+1), indexlai) |
---|
2196 | CALL histwrite_p(hist_id, 'assimi', kjit, assimi, kjpindex*(nlai+1), indexlai) |
---|
2197 | CALL histwrite_p(hist_id, 'Rd', kjit, Rd, kjpindex*(nlai+1), indexlai) |
---|
2198 | ENDIF |
---|
2199 | !! Calculated intercellular CO2 over nlai needed for the chemistry module |
---|
2200 | cim(:,jv)=0. |
---|
2201 | laisum(:)=0 |
---|
2202 | DO jl=1,nlai |
---|
2203 | WHERE (laitab(jl) .LE. lai(:,jv) ) |
---|
2204 | cim(:,jv)= cim(:,jv)+leaf_ci(:,jv,jl)*(laitab(jl+1)-laitab(jl)) |
---|
2205 | laisum(:)=laisum(:)+ (laitab(jl+1)-laitab(jl)) |
---|
2206 | ENDWHERE |
---|
2207 | ENDDO |
---|
2208 | WHERE (laisum(:)>0) |
---|
2209 | cim(:,jv)= cim(:,jv)/laisum(:) |
---|
2210 | ENDWHERE |
---|
2211 | |
---|
2212 | |
---|
2213 | ! |
---|
2214 | !! 2.5 Calculate resistances |
---|
2215 | ! |
---|
2216 | IF (nia .GT. 0) THEN |
---|
2217 | ! |
---|
2218 | DO inia=1,nia |
---|
2219 | ! |
---|
2220 | iainia=index_assi(inia) |
---|
2221 | |
---|
2222 | !! Mean stomatal conductance for CO2 (mol m-2 s-1) |
---|
2223 | gsmean(iainia,jv) = gstot(iainia) |
---|
2224 | ! |
---|
2225 | ! cimean is the "mean ci" calculated in such a way that assimilation |
---|
2226 | ! calculated in enerbil is equivalent to assimtot |
---|
2227 | ! |
---|
2228 | IF ( ABS(gsmean(iainia,jv)-g0var(iainia)*laisum(iainia)) .GT. min_sechiba) THEN |
---|
2229 | cimean(iainia,jv) = (fvpd(iainia)*(assimtot(iainia)+Rdtot(iainia))) /& |
---|
2230 | (gsmean(iainia,jv)-g0var(iainia)*laisum(iainia)) + gamma_star(iainia) |
---|
2231 | ELSE |
---|
2232 | cimean(iainia,jv) = gamma_star(iainia) |
---|
2233 | ENDIF |
---|
2234 | |
---|
2235 | ! conversion from umol m-2 (PFT) s-1 to gC m-2 (mesh area) tstep-1 |
---|
2236 | gpp(iainia,jv) = assimtot(iainia)*12e-6*veget_max(iainia,jv)*dt_sechiba |
---|
2237 | |
---|
2238 | ! |
---|
2239 | ! conversion from mol/m^2/s to m/s |
---|
2240 | ! |
---|
2241 | ! As in Pearcy, Schulze and Zimmermann |
---|
2242 | ! Measurement of transpiration and leaf conductance |
---|
2243 | ! Chapter 8 of Plant Physiological Ecology |
---|
2244 | ! Field methods and instrumentation, 1991 |
---|
2245 | ! Editors: |
---|
2246 | ! |
---|
2247 | ! Robert W. Pearcy, |
---|
2248 | ! James R. Ehleringer, |
---|
2249 | ! Harold A. Mooney, |
---|
2250 | ! Philip W. Rundel |
---|
2251 | ! |
---|
2252 | ! ISBN: 978-0-412-40730-7 (Print) 978-94-010-9013-1 (Online) |
---|
2253 | |
---|
2254 | gstot(iainia) = mol_to_m_1 *(temp_air(iainia)/tp_00)*& |
---|
2255 | (pb_std/pb(iainia))*gstot(iainia)*ratio_H2O_to_CO2 |
---|
2256 | gstop(iainia) = mol_to_m_1 * (temp_air(iainia)/tp_00)*& |
---|
2257 | (pb_std/pb(iainia))*leaf_gs_top(iainia)*ratio_H2O_to_CO2*& |
---|
2258 | laitab(ilai(iainia)+1) |
---|
2259 | ! |
---|
2260 | rveget(iainia,jv) = un/gstop(iainia) |
---|
2261 | |
---|
2262 | ! |
---|
2263 | ! |
---|
2264 | ! rstruct is the difference between rtot (=1./gstot) and rveget |
---|
2265 | ! |
---|
2266 | ! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif |
---|
2267 | !rstruct(iainia,jv) = un/gstot(iainia) - & |
---|
2268 | ! rveget(iainia,jv) |
---|
2269 | rstruct(iainia,jv) = MAX( un/gstot(iainia) - & |
---|
2270 | rveget(iainia,jv), min_sechiba) |
---|
2271 | ! |
---|
2272 | ! |
---|
2273 | !! wind is a global variable of the diffuco module. |
---|
2274 | speed = MAX(min_wind, wind(iainia)) |
---|
2275 | ! |
---|
2276 | ! beta for transpiration |
---|
2277 | ! |
---|
2278 | ! Corrections Nathalie - 28 March 2006 - on advices of Fred Hourdin |
---|
2279 | !! Introduction of a potentiometer rveg_pft to settle the rveg+rstruct sum problem in the coupled mode. |
---|
2280 | !! rveg_pft=1 in the offline mode. rveg_pft is a global variable declared in the diffuco module. |
---|
2281 | !vbeta3(iainia,jv) = veget_max(iainia,jv) * & |
---|
2282 | ! (un - zqsvegrap(iainia)) * & |
---|
2283 | ! (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + & |
---|
2284 | ! rstruct(iainia,jv)))) |
---|
2285 | !! Global resistance of the canopy to evaporation |
---|
2286 | cresist=(un / (un + speed * q_cdrag(iainia) * & |
---|
2287 | veget(iainia,jv)/veget_max(iainia,jv) * & |
---|
2288 | (rveg_pft(jv)*(rveget(iainia,jv) + rstruct(iainia,jv))))) |
---|
2289 | |
---|
2290 | IF ( humrel(iainia,jv) >= min_sechiba ) THEN |
---|
2291 | vbeta3(iainia,jv) = veget(iainia,jv) * & |
---|
2292 | (un - zqsvegrap(iainia)) * cresist + & |
---|
2293 | MIN( vbeta23(iainia,jv), veget(iainia,jv) * & |
---|
2294 | zqsvegrap(iainia) * cresist ) |
---|
2295 | ELSE |
---|
2296 | ! Because of a minimum conductance g0, vbeta3 cannot be zero even if humrel=0 |
---|
2297 | ! in the above equation. |
---|
2298 | ! Here, we force transpiration to be zero when the soil cannot deliver it |
---|
2299 | vbeta3(iainia,jv) = zero |
---|
2300 | END IF |
---|
2301 | |
---|
2302 | ! vbeta3pot for computation of potential transpiration (needed for irrigation) |
---|
2303 | vbeta3pot(iainia,jv) = MAX(zero, veget(iainia,jv) * cresist) |
---|
2304 | ! |
---|
2305 | ! |
---|
2306 | ENDDO |
---|
2307 | ! |
---|
2308 | ENDIF |
---|
2309 | ! |
---|
2310 | END DO ! loop over vegetation types |
---|
2311 | ! |
---|
2312 | |
---|
2313 | |
---|
2314 | IF (printlev>=3) WRITE (numout,*) ' diffuco_trans_co2 done ' |
---|
2315 | |
---|
2316 | |
---|
2317 | END SUBROUTINE diffuco_trans_co2 |
---|
2318 | |
---|
2319 | |
---|
2320 | !! ================================================================================================================================ |
---|
2321 | !! SUBROUTINE : diffuco_comb |
---|
2322 | !! |
---|
2323 | !>\BRIEF This routine combines the previous partial beta |
---|
2324 | !! coefficients and calculates the total alpha and complete beta coefficients. |
---|
2325 | !! |
---|
2326 | !! DESCRIPTION : Those integrated coefficients are used to calculate (in enerbil.f90) the total evapotranspiration |
---|
2327 | !! from the grid-cell. \n |
---|
2328 | !! |
---|
2329 | !! In the case that air is more humid than surface, dew deposition can occur (negative latent heat flux). |
---|
2330 | !! In this instance, for temperature above zero, all of the beta coefficients are set to 0, except for |
---|
2331 | !! interception (vbeta2) and bare soil (vbeta4 with zero soil resistance). The amount of water that is |
---|
2332 | !! intercepted by leaves is calculated based on the value of LAI of the surface. In the case of freezing |
---|
2333 | !! temperatures, water is added to the snow reservoir, and so vbeta4 and vbeta2 are set to 0, and the |
---|
2334 | !! total vbeta is set to 1.\n |
---|
2335 | !! |
---|
2336 | !! \latexonly |
---|
2337 | !! \input{diffucocomb1.tex} |
---|
2338 | !! \endlatexonly |
---|
2339 | !! |
---|
2340 | !! The beta and alpha coefficients are initially set to 1. |
---|
2341 | !! \latexonly |
---|
2342 | !! \input{diffucocomb2.tex} |
---|
2343 | !! \endlatexonly |
---|
2344 | !! |
---|
2345 | !! If snow is lower than the critical value: |
---|
2346 | !! \latexonly |
---|
2347 | !! \input{diffucocomb3.tex} |
---|
2348 | !! \endlatexonly |
---|
2349 | !! If in the presence of dew: |
---|
2350 | !! \latexonly |
---|
2351 | !! \input{diffucocomb4.tex} |
---|
2352 | !! \endlatexonly |
---|
2353 | !! |
---|
2354 | !! Determine where the water goes (soil, vegetation, or snow) |
---|
2355 | !! when air moisture exceeds saturation. |
---|
2356 | !! \latexonly |
---|
2357 | !! \input{diffucocomb5.tex} |
---|
2358 | !! \endlatexonly |
---|
2359 | !! |
---|
2360 | !! If it is not freezing dew is put into the interception reservoir and onto the bare soil. If it is freezing, |
---|
2361 | !! water is put into the snow reservoir. |
---|
2362 | !! Now modify vbetas where necessary: for soil and snow |
---|
2363 | !! \latexonly |
---|
2364 | !! \input{diffucocomb6.tex} |
---|
2365 | !! \endlatexonly |
---|
2366 | !! |
---|
2367 | !! and for vegetation |
---|
2368 | !! \latexonly |
---|
2369 | !! \input{diffucocomb7.tex} |
---|
2370 | !! \endlatexonly |
---|
2371 | !! |
---|
2372 | !! Then compute part of dew that can be intercepted by leafs. |
---|
2373 | !! |
---|
2374 | !! There will be no transpiration when air moisture is too high, under any circumstance |
---|
2375 | !! \latexonly |
---|
2376 | !! \input{diffucocomb8.tex} |
---|
2377 | !! \endlatexonly |
---|
2378 | !! |
---|
2379 | !! There will also be no interception loss on bare soil, under any circumstance. |
---|
2380 | !! \latexonly |
---|
2381 | !! \input{diffucocomb9.tex} |
---|
2382 | !! \endlatexonly |
---|
2383 | !! |
---|
2384 | !! The flowchart details the 'decision tree' which underlies the module. |
---|
2385 | !! |
---|
2386 | !! RECENT CHANGE(S): None |
---|
2387 | !! |
---|
2388 | !! MAIN OUTPUT VARIABLE(S): vbeta1, vbeta4, humrel, vbeta2, vbeta3, vbeta |
---|
2389 | !! |
---|
2390 | !! REFERENCE(S) : |
---|
2391 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
2392 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
2393 | !! Circulation Model. Journal of Climate, 6, pp.248-273 |
---|
2394 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation |
---|
2395 | !! sur le cycle de l'eau, PhD Thesis, available from: |
---|
2396 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf |
---|
2397 | !! |
---|
2398 | !! FLOWCHART : |
---|
2399 | !! \latexonly |
---|
2400 | !! \includegraphics[scale=0.25]{diffuco_comb_flowchart.png} |
---|
2401 | !! \endlatexonly |
---|
2402 | !! \n |
---|
2403 | !_ ================================================================================================================================ |
---|
2404 | |
---|
2405 | SUBROUTINE diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, & |
---|
2406 | & snow, veget, lai, tot_bare_soil, vbeta1, vbeta2, vbeta3 , vbeta4, & |
---|
2407 | & evap_bare_lim, evap_bare_lim_ns, veget_max, vbeta, qsintmax) |
---|
2408 | |
---|
2409 | ! Ajout qsintmax dans les arguments de la routine Nathalie / le 13-03-2006 |
---|
2410 | |
---|
2411 | !! 0. Variable and parameter declaration |
---|
2412 | |
---|
2413 | !! 0.1 Input variables |
---|
2414 | |
---|
2415 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
2416 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3}) |
---|
2417 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1}) |
---|
2418 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Nortward Lowest level wind speed (m s^{-1}) |
---|
2419 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
2420 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa) |
---|
2421 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1}) |
---|
2422 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature (K) |
---|
2423 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lower air temperature (K) |
---|
2424 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg) |
---|
2425 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (fraction) |
---|
2426 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2}) |
---|
2427 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2}) |
---|
2428 | REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil!! Total evaporating bare soil fraction |
---|
2429 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty) |
---|
2430 | |
---|
2431 | !! 0.2 Output variables |
---|
2432 | |
---|
2433 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient (-) |
---|
2434 | |
---|
2435 | !! 0.3 Modified variables |
---|
2436 | |
---|
2437 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta1 !! Beta for sublimation (-) |
---|
2438 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta4 !! Beta for Bare soil evaporation (-) |
---|
2439 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Soil moisture stress (within range 0 to 1) |
---|
2440 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta2 !! Beta for interception loss (-) |
---|
2441 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3 !! Beta for Transpiration (-) |
---|
2442 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evap_bare_lim !! limiting factor for bare soil evaporation |
---|
2443 | !! when the 11-layer hydrology is used (-) |
---|
2444 | REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout):: evap_bare_lim_ns !! limiting factor for bare soil evaporation |
---|
2445 | !! when the 11-layer hydrology is used (-) |
---|
2446 | !! 0.4 Local variables |
---|
2447 | |
---|
2448 | INTEGER(i_std) :: ji, jv |
---|
2449 | REAL(r_std) :: zevtest, zsoil_moist, zrapp |
---|
2450 | REAL(r_std), DIMENSION(kjpindex) :: qsatt |
---|
2451 | LOGICAL, DIMENSION(kjpindex) :: toveg, tosnow |
---|
2452 | REAL(r_std) :: coeff_dew_veg |
---|
2453 | REAL(r_std), DIMENSION(kjpindex) :: vegtot |
---|
2454 | !_ ================================================================================================================================ |
---|
2455 | |
---|
2456 | !! 1 If we are in presence of dew |
---|
2457 | |
---|
2458 | CALL qsatcalc (kjpindex, temp_sol, pb, qsatt) |
---|
2459 | |
---|
2460 | |
---|
2461 | !! 1.1 Determine where the water goes |
---|
2462 | !! Determine where the water goes (soil, vegetation, or snow) |
---|
2463 | !! when air moisture exceeds saturation. |
---|
2464 | !! \latexonly |
---|
2465 | !! \input{diffucocomb5.tex} |
---|
2466 | !! \endlatexonly |
---|
2467 | toveg(:) = .FALSE. |
---|
2468 | tosnow(:) = .FALSE. |
---|
2469 | DO ji = 1, kjpindex |
---|
2470 | IF ( qsatt(ji) .LT. qair(ji) ) THEN |
---|
2471 | IF (temp_air(ji) .GT. tp_00) THEN |
---|
2472 | !! If it is not freezing dew is put into the |
---|
2473 | !! interception reservoir and onto the bare soil. |
---|
2474 | toveg(ji) = .TRUE. |
---|
2475 | ELSE |
---|
2476 | !! If it is freezing water is put into the |
---|
2477 | !! snow reservoir. |
---|
2478 | tosnow(ji) = .TRUE. |
---|
2479 | ENDIF |
---|
2480 | ENDIF |
---|
2481 | END DO |
---|
2482 | |
---|
2483 | !! 1.2 Now modify vbetas where necessary. |
---|
2484 | |
---|
2485 | !! 1.2.1 Soil and snow |
---|
2486 | !! \latexonly |
---|
2487 | !! \input{diffucocomb6.tex} |
---|
2488 | !! \endlatexonly |
---|
2489 | |
---|
2490 | ! We need to keep consistency between evap_bare_lim, evap_bare_lim_ns and vbeta4 (thus vevapnu) |
---|
2491 | ! or we have a water conservation issue in hydrol_split_soil |
---|
2492 | |
---|
2493 | DO ji = 1, kjpindex |
---|
2494 | |
---|
2495 | IF ( toveg(ji) ) THEN |
---|
2496 | |
---|
2497 | vbeta1(ji) = zero |
---|
2498 | vegtot(ji) = SUM(veget_max(ji,:)) |
---|
2499 | |
---|
2500 | IF ( (tot_bare_soil(ji) .GT. min_sechiba) .AND. (vegtot(ji).GT. min_sechiba) ) THEN |
---|
2501 | |
---|
2502 | vbeta4(ji) = tot_bare_soil(ji) |
---|
2503 | |
---|
2504 | evap_bare_lim_ns(ji,:) = tot_bare_soil(ji)/vegtot(ji) |
---|
2505 | |
---|
2506 | evap_bare_lim(ji) = vbeta4(ji) |
---|
2507 | ! consistent with evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)*soiltile(ji,:)*vegtot(ji)) |
---|
2508 | ! as SUM(soiltile(ji,:)) = 1 |
---|
2509 | |
---|
2510 | ELSE |
---|
2511 | vbeta4(ji) = zero |
---|
2512 | evap_bare_lim_ns(ji,:) = zero |
---|
2513 | evap_bare_lim(ji) = zero |
---|
2514 | ENDIF |
---|
2515 | ENDIF |
---|
2516 | |
---|
2517 | IF ( tosnow(ji) ) THEN |
---|
2518 | vbeta1(ji) = un |
---|
2519 | vbeta4(ji) = zero |
---|
2520 | evap_bare_lim_ns(ji,:) = zero |
---|
2521 | evap_bare_lim(ji) = zero |
---|
2522 | ENDIF |
---|
2523 | |
---|
2524 | ENDDO |
---|
2525 | |
---|
2526 | !! 1.2.2 Vegetation and interception loss |
---|
2527 | !! \latexonly |
---|
2528 | !! \input{diffucocomb7.tex} |
---|
2529 | !! \endlatexonly |
---|
2530 | DO jv = 1, nvm |
---|
2531 | |
---|
2532 | DO ji = 1, kjpindex |
---|
2533 | |
---|
2534 | IF ( toveg(ji) ) THEN |
---|
2535 | IF (qsintmax(ji,jv) .GT. min_sechiba) THEN |
---|
2536 | |
---|
2537 | ! Compute part of dew that can be intercepted by leafs. |
---|
2538 | IF ( lai(ji,jv) .GT. min_sechiba) THEN |
---|
2539 | IF (lai(ji,jv) .GT. 1.5) THEN |
---|
2540 | coeff_dew_veg= & |
---|
2541 | & dew_veg_poly_coeff(6)*lai(ji,jv)**5 & |
---|
2542 | & - dew_veg_poly_coeff(5)*lai(ji,jv)**4 & |
---|
2543 | & + dew_veg_poly_coeff(4)*lai(ji,jv)**3 & |
---|
2544 | & - dew_veg_poly_coeff(3)*lai(ji,jv)**2 & |
---|
2545 | & + dew_veg_poly_coeff(2)*lai(ji,jv) & |
---|
2546 | & + dew_veg_poly_coeff(1) |
---|
2547 | ELSE |
---|
2548 | coeff_dew_veg=un |
---|
2549 | ENDIF |
---|
2550 | ELSE |
---|
2551 | coeff_dew_veg=zero |
---|
2552 | ENDIF |
---|
2553 | IF (jv .EQ. 1) THEN |
---|
2554 | ! This line may not work with CWRR when frac_bare is distributed among three soiltiles |
---|
2555 | ! Fortunately, qsintmax(ji,1)=0 (LAI=0 in PFT1) so we never pass here |
---|
2556 | vbeta2(ji,jv) = coeff_dew_veg*tot_bare_soil(ji) |
---|
2557 | ELSE |
---|
2558 | vbeta2(ji,jv) = coeff_dew_veg*veget(ji,jv) |
---|
2559 | ENDIF |
---|
2560 | ELSE |
---|
2561 | vbeta2(ji,jv) = zero ! if qsintmax=0, vbeta2=0 |
---|
2562 | ENDIF |
---|
2563 | ENDIF |
---|
2564 | IF ( tosnow(ji) ) vbeta2(ji,jv) = zero |
---|
2565 | |
---|
2566 | ENDDO |
---|
2567 | |
---|
2568 | ENDDO |
---|
2569 | |
---|
2570 | !! 1.2.3 Vegetation and transpiration |
---|
2571 | !! There will be no transpiration when air moisture is too high, under any circumstance |
---|
2572 | !! \latexonly |
---|
2573 | !! \input{diffucocomb8.tex} |
---|
2574 | !! \endlatexonly |
---|
2575 | DO jv = 1, nvm |
---|
2576 | DO ji = 1, kjpindex |
---|
2577 | IF ( qsatt(ji) .LT. qair(ji) ) THEN |
---|
2578 | vbeta3(ji,jv) = zero |
---|
2579 | humrel(ji,jv) = zero |
---|
2580 | ENDIF |
---|
2581 | ENDDO |
---|
2582 | ENDDO |
---|
2583 | |
---|
2584 | |
---|
2585 | !! 1.2.4 Overrules 1.2.2 |
---|
2586 | !! There will also be no interception loss on bare soil, under any circumstance. |
---|
2587 | !! \latexonly |
---|
2588 | !! \input{diffucocomb9.tex} |
---|
2589 | !! \endlatexonly |
---|
2590 | DO ji = 1, kjpindex |
---|
2591 | IF ( qsatt(ji) .LT. qair(ji) ) THEN |
---|
2592 | vbeta2(ji,1) = zero |
---|
2593 | ENDIF |
---|
2594 | ENDDO |
---|
2595 | |
---|
2596 | !! 2 Now calculate vbeta in all cases (the equality needs to hold for enerbil to be consistent) |
---|
2597 | |
---|
2598 | DO ji = 1, kjpindex |
---|
2599 | vbeta(ji) = vbeta4(ji) + SUM(vbeta2(ji,:)) + SUM(vbeta3(ji,:)) |
---|
2600 | |
---|
2601 | IF (vbeta(ji) .LT. min_sechiba) THEN |
---|
2602 | vbeta(ji) = zero |
---|
2603 | vbeta4(ji) = zero |
---|
2604 | vbeta2(ji,:)= zero |
---|
2605 | vbeta3(ji,:)= zero |
---|
2606 | evap_bare_lim_ns(ji,:) = zero |
---|
2607 | evap_bare_lim(ji) = zero |
---|
2608 | END IF |
---|
2609 | ENDDO |
---|
2610 | |
---|
2611 | CALL xios_orchidee_send_field("evap_bare_lim",evap_bare_lim) |
---|
2612 | CALL xios_orchidee_send_field("evap_bare_lim_ns",evap_bare_lim_ns) |
---|
2613 | |
---|
2614 | IF (printlev>=3) WRITE (numout,*) ' diffuco_comb done ' |
---|
2615 | |
---|
2616 | END SUBROUTINE diffuco_comb |
---|
2617 | |
---|
2618 | |
---|
2619 | !! ================================================================================================================================ |
---|
2620 | !! SUBROUTINE : diffuco_raerod |
---|
2621 | !! |
---|
2622 | !>\BRIEF Computes the aerodynamic resistance, for cases in which the |
---|
2623 | !! surface drag coefficient is provided by the coupled atmospheric model LMDZ and when the flag |
---|
2624 | !! 'ldq_cdrag_from_gcm' is set to TRUE |
---|
2625 | !! |
---|
2626 | !! DESCRIPTION : Simply computes the aerodynamic resistance, for cases in which the |
---|
2627 | !! surface drag coefficient is provided by the coupled atmospheric model LMDZ. If the surface drag coefficient |
---|
2628 | !! is not provided by the LMDZ or signalled by the flag 'ldq_cdrag_from_gcm' set to FALSE, then the subroutine |
---|
2629 | !! diffuco_aero is called instead of this one. |
---|
2630 | !! |
---|
2631 | !! Calculation of the aerodynamic resistance, for diganostic purposes. First calculate wind speed: |
---|
2632 | !! \latexonly |
---|
2633 | !! \input{diffucoaerod1.tex} |
---|
2634 | !! \endlatexonly |
---|
2635 | !! |
---|
2636 | !! next calculate ::raero |
---|
2637 | !! \latexonly |
---|
2638 | !! \input{diffucoaerod2.tex} |
---|
2639 | !! \endlatexonly |
---|
2640 | !! |
---|
2641 | !! RECENT CHANGE(S): None |
---|
2642 | !! |
---|
2643 | !! MAIN OUTPUT VARIABLE(S): ::raero |
---|
2644 | !! |
---|
2645 | !! REFERENCE(S) : |
---|
2646 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
2647 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
2648 | !! Circulation Model. Journal of Climate, 6, pp.248-273 |
---|
2649 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influence de l'irrigation |
---|
2650 | !! sur le cycle de l'eau, PhD Thesis, available from: |
---|
2651 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf |
---|
2652 | !! |
---|
2653 | !! FLOWCHART : None |
---|
2654 | !! \n |
---|
2655 | !_ ================================================================================================================================ |
---|
2656 | |
---|
2657 | SUBROUTINE diffuco_raerod (kjpindex, u, v, q_cdrag, raero) |
---|
2658 | |
---|
2659 | IMPLICIT NONE |
---|
2660 | |
---|
2661 | !! 0. Variable and parameter declaration |
---|
2662 | |
---|
2663 | !! 0.1 Input variables |
---|
2664 | |
---|
2665 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
2666 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind velocity (m s^{-1}) |
---|
2667 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind velocity (m s^{-1}) |
---|
2668 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
2669 | |
---|
2670 | !! 0.2 Output variables |
---|
2671 | |
---|
2672 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: raero !! Aerodynamic resistance (s m^{-1}) |
---|
2673 | |
---|
2674 | !! 0.3 Modified variables |
---|
2675 | |
---|
2676 | !! 0.4 Local variables |
---|
2677 | |
---|
2678 | INTEGER(i_std) :: ji !! (-) |
---|
2679 | REAL(r_std) :: speed !! (m s^{-1}) |
---|
2680 | !_ ================================================================================================================================ |
---|
2681 | |
---|
2682 | !! 1. Simple calculation of the aerodynamic resistance, for diganostic purposes. |
---|
2683 | |
---|
2684 | DO ji=1,kjpindex |
---|
2685 | |
---|
2686 | !! \latexonly |
---|
2687 | !! \input{diffucoaerod1.tex} |
---|
2688 | !! \endlatexonly |
---|
2689 | speed = MAX(min_wind, wind(ji)) |
---|
2690 | |
---|
2691 | !! \latexonly |
---|
2692 | !! \input{diffucoaerod2.tex} |
---|
2693 | !! \endlatexonly |
---|
2694 | raero(ji) = un / (q_cdrag(ji)*speed) |
---|
2695 | |
---|
2696 | ENDDO |
---|
2697 | |
---|
2698 | END SUBROUTINE diffuco_raerod |
---|
2699 | |
---|
2700 | |
---|
2701 | FUNCTION Arrhenius (kjpindex,temp,ref_temp,energy_act) RESULT ( val_arrhenius ) |
---|
2702 | !! 0.1 Input variables |
---|
2703 | |
---|
2704 | INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-) |
---|
2705 | REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K) |
---|
2706 | REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K) |
---|
2707 | REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1) |
---|
2708 | |
---|
2709 | !! 0.2 Result |
---|
2710 | |
---|
2711 | REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on |
---|
2712 | !! a Arrhenius function (-) |
---|
2713 | |
---|
2714 | val_arrhenius(:)=EXP(((temp(:)-ref_temp)*energy_act)/(ref_temp*RR*(temp(:)))) |
---|
2715 | END FUNCTION Arrhenius |
---|
2716 | |
---|
2717 | FUNCTION Arrhenius_modified_1d (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius ) |
---|
2718 | !! 0.1 Input variables |
---|
2719 | |
---|
2720 | INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-) |
---|
2721 | REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K) |
---|
2722 | REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K) |
---|
2723 | REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1) |
---|
2724 | REAL(r_std),INTENT(in) :: energy_deact !! Deactivation Energy (J mol-1) |
---|
2725 | REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: entropy !! Entropy term (J K-1 mol-1) |
---|
2726 | |
---|
2727 | !! 0.2 Result |
---|
2728 | |
---|
2729 | REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on |
---|
2730 | !! a Arrhenius function (-) |
---|
2731 | |
---|
2732 | val_arrhenius(:)=EXP(((temp(:)-ref_temp)*energy_act)/(ref_temp*RR*(temp(:)))) & |
---|
2733 | * (1. + EXP( (ref_temp * entropy(:) - energy_deact) / (ref_temp * RR ))) & |
---|
2734 | / (1. + EXP( (temp(:) * entropy(:) - energy_deact) / ( RR*temp(:)))) |
---|
2735 | |
---|
2736 | END FUNCTION Arrhenius_modified_1d |
---|
2737 | |
---|
2738 | FUNCTION Arrhenius_modified_0d (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius ) |
---|
2739 | !! 0.1 Input variables |
---|
2740 | |
---|
2741 | INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-) |
---|
2742 | REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K) |
---|
2743 | REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K) |
---|
2744 | REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1) |
---|
2745 | REAL(r_std),INTENT(in) :: energy_deact !! Deactivation Energy (J mol-1) |
---|
2746 | REAL(r_std),INTENT(in) :: entropy !! Entropy term (J K-1 mol-1) |
---|
2747 | |
---|
2748 | !! 0.2 Result |
---|
2749 | |
---|
2750 | REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on |
---|
2751 | !! a Arrhenius function (-) |
---|
2752 | |
---|
2753 | val_arrhenius(:)=EXP(((temp(:)-ref_temp)*energy_act)/(ref_temp*RR*(temp(:)))) & |
---|
2754 | * (1. + EXP( (ref_temp * entropy - energy_deact) / (ref_temp * RR ))) & |
---|
2755 | / (1. + EXP( (temp(:) * entropy - energy_deact) / ( RR*temp(:)))) |
---|
2756 | |
---|
2757 | END FUNCTION Arrhenius_modified_0d |
---|
2758 | |
---|
2759 | |
---|
2760 | END MODULE diffuco |
---|