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