MODULE FARQUHAR_MOD USE PARKIND1, ONLY : JPIM, JPRB, JPRM USE YOMHOOK, ONLY : LHOOK, DR_HOOK USE YOS_AGS, ONLY : TAGS USE YOS_AGF, ONLY : TAGF USE SUFARQUHAR_MOD IMPLICIT NONE INTERFACE Arrhenius_modified MODULE PROCEDURE Arrhenius_modified_0d, Arrhenius_modified_1d END INTERFACE REAL(KIND=JPRB), ALLOCATABLE, SAVE :: temp_growth(:) ! Growth temperature (°C) - Is equal to t2m_month LOGICAL, SAVE :: FIRST_CALL=.TRUE. INTEGER(KIND=JPIM), SAVE :: KSTEP_SAV=-1 CONTAINS !! ============================================================================================================================== !! SUBROUTINE : FARQUHAR !! !! This subroutine computes carbon assimilation and stomatal !! conductance, following respectively Farqhuar et al. (1980) and Ball et al. (1987). !! !! DESCRIPTION: !! *** General: !! The equations are different depending on the photosynthesis mode (C3 versus C4). !! Assimilation and conductance are computed over 20 levels of LAI and then !! integrated at the canopy level. !! This routine also computes partial beta coefficient: transpiration for each !! type of vegetation. !! There is a main loop on the PFTs, then inner loops on the points where !! assimilation has to be calculated. !! This subroutine is called at each sechiba time step by sechiba_main. !! *** Details: !! - Integration at the canopy level !! - Light's extinction !! The available light follows a simple Beer extinction law. !! The extinction coefficients (ext_coef) are PFT-dependant constants. !! - Estimation of relative humidity of air (for calculation of the stomatal conductance) !! - Calculation of the water limitation factor !! - Calculation of temperature dependent parameters for C4 plants !! - Calculation of temperature dependent parameters for C3 plants !! - Vmax scaling !! Vmax is scaled into the canopy due to reduction of nitrogen !! (Johnson and Thornley,1984). !! - Assimilation for C4 plants (Collatz et al., 1992) !! - Assimilation for C3 plants (Farqhuar et al., 1980) !! - Estimation of the stomatal conductance (Ball et al., 1987) !! !! !! REFERENCE(S) : !! - Ball, J., T. Woodrow, and J. Berry (1987), A model predicting stomatal !! conductance and its contribution to the control of photosynthesis under !! different environmental conditions, Prog. Photosynthesis, 4, 221– 224. !! - Collatz, G., M. Ribas-Carbo, and J. Berry (1992), Coupled photosynthesis !! stomatal conductance model for leaves of C4 plants, Aust. J. Plant Physiol., !! 19, 519–538. !! - Farquhar, G., S. von Caemmener, and J. Berry (1980), A biochemical model of !! photosynthesis CO2 fixation in leaves of C3 species, Planta, 149, 78–90. !! - Johnson, I. R., and J. Thornley (1984), A model of instantaneous and daily !! canopy photosynthesis, J Theor. Biol., 107, 531 545 !! - McMurtrie, R.E., Rook, D.A. and Kelliher, F.M., 1990. Modelling the yield of Pinus radiata on a !! site limited by water and nitrogen. For. Ecol. Manage., 30: 381-413 !! - Bounoua, L., Hall, F. G., Sellers, P. J., Kumar, A., Collatz, G. J., Tucker, C. J., and Imhoff, M. L. (2010), Quantifying the !! negative feedback of vegetation to greenhouse warming: A modeling approach, Geophysical Research Letters, 37, Artn L23701, !! Doi 10.1029/2010gl045338 !! - Bounoua, L., Collatz, G. J., Sellers, P. J., Randall, D. A., Dazlich, D. A., Los, S. O., Berry, J. A., Fung, I., !! Tucker, C. J., Field, C. B., and Jensen, T. G. (1999), Interactions between vegetation and climate: Radiative and physiological !! effects of doubled atmospheric co2, Journal of Climate, 12, 309-324, Doi 10.1175/1520-0442(1999)012<0309:Ibvacr>2.0.Co;2 !! - Sellers, P. J., Bounoua, L., Collatz, G. J., Randall, D. A., Dazlich, D. A., Los, S. O., Berry, J. A., Fung, I., !! Tucker, C. J., Field, C. B., and Jensen, T. G. (1996), Comparison of radiative and physiological effects of doubled atmospheric !! co2 on climate, Science, 271, 1402-1406, DOI 10.1126/science.271.5254.1402 !! - Lewis, J. D., Ward, J. K., and Tissue, D. T. (2010), Phosphorus supply drives nonlinear responses of cottonwood !! (populus deltoides) to increases in co2 concentration from glacial to future concentrations, New Phytologist, 187, 438-448, !! DOI 10.1111/j.1469-8137.2010.03307.x !! - Kattge, J., Knorr, W., Raddatz, T., and Wirth, C. (2009), Quantifying photosynthetic capacity and its relationship to leaf !! nitrogen content for global-scale terrestrial biosphere models, Global Change Biology, 15, 976-991, !! DOI 10.1111/j.1365-2486.2008.01744.x !! !_ ================================================================================================================================ !Initial ORCHIDEE SUBROUTINE !SUBROUTINE diffuco_trans_co2 (kjpindex, swdown, pb, qsurf, qair, temp_air, temp_growth, rau, u, v, q_cdrag, humrel, & ! assim_param, Ca, & ! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, vbeta3pot, rveget, rstruct, & ! cimean, gsmean, gpp, & ! co2_to_bm, vbeta23, index, kjit, cim) !CALL in cotworestress_mod.F90 !CALL FARQUHAR(KIDIA, KFDIA, KLON, LDLAND, KVTYPE, YDAGF, KSTEP, PTSTEP, PSRFD, PAPHM1/100., PQS, PQM1, PTM1, PTSKM1M, PF2, & ! & RCO2/44.*29.E6, PLAI, ZXTGS, PAGF) SUBROUTINE FARQUHAR (KIDIA, KFDIA, KLON, LDLAND, KVTYPE, YDAGF, KSTEP, dt_sechiba, swdown, pb, qsurf, qair, temp_air, skin_temp, humrel, & Ca, lai, gsmean, gpp) ! INTERFACE ! --------- ! ! INPUT PARAMETERS ! ---------------- ! KIDIA Begin point in arrays (from surfexcdriver_ctl_mod.F90) ! KFDIA End point in arrays (from surfexcdriver_ctl_mod.F90) ! KLON Length of arrays (from surfexcdriver_ctl_mod.F90) ! KVTYPE VEGETATION TYPE CORRESPONDING TO TILE (from cotworestress_mod.F90) ! LDLAND (from cotworestress_mod.F90) ! KSTEP Time step index (from surfexcdriver_ctl_mod.F90) ! PTSTEP Timestep (from surfexcdriver_ctl_mod.F90) ! PTM1 TEMPERATURE AT T-1 K (from cotworestress_mod.F90) ! PTSKM1M SURFACE TEMPERATURE K (from cotworestress_mod.F90) ! PQM1 SPECIFIC HUMIDITY AT T-1 KG/KG (from cotworestress_mod.F90) ! PAPHM1 PRESSURE AT T-1 PA (from cotworestress_mod.F90) !!!@FM change unit to hPa ! PLAI LEAF AREA INDEX M2/M2 (from cotworestress_mod.F90) ! PSRFD DOWNWARD SHORT WAVE RADIATION FLUX AT SURFACE W/M2 (from cotworestress_mod.F90) ! PQS SATURATION Q AT SURFACE KG/KG (from cotworestress_mod.F90) ! PF2 SOIL MOISTURE STRESS FUNCTION - (from cotworestress_mod.F90) ! RCO2 atmospheric CO2 concentration (kgCO2 kgAir-1) (from yos_ags.F90) !!!@FM change unit to ppm ! ! OUTPUT PARAMETERS ! ----------------- ! PAGF GROSS ASSIMILATION OF CO2 KG_CO2 KG_AIR-1 M S-1 ! PGSF LEAF CONDUCTANCE TO H20 (CUTICULAR AND STOMATAL) M/S ! ! ! MATCHING TABLE WITH ORCHIDEE PFTS ! --------------------------------- ! (1) ! Crops, Mixed Farming => PFT13 (C4 crops) ! (2) ! Short Grass => PFT10 (C3 grass) ! (3) ! Evergreen Needleleaf Trees => PFT7 (Boreal Needleleaf Evergreen) ! (4) ! Deciduous Needleleaf Trees => PFT9 (Boreal Needleleaf Deciduous) ! (5) ! Deciduous Broadleaf Trees => PFT6 (Temperate Broadleaf Summergreen) ! (6) ! Evergreen Broadleaf Trees => PFT5 (Temperate Broadleaf Evergreen) ! (7) ! Tall Grass => PFT11 (C4 grass) ! (8) ! Desert => none ! (9) ! Tundra => PFT10 (C3 grass) ! (10) ! Irrigated Crops => PFT12 (C3 crops) ! (11) ! Semidesert => PFT11 (C4 grass) ! (12) ! Ice Caps and Glaciers => none ! (13) ! Bogs and Marshes => PFT10 (C3 grass) ! (14) ! Inland Water => none ! (15) ! Ocean => none ! (16) ! Evergreen Shrubs => PFT11 (C4 grass) ! (17) ! Deciduous Shrubs => PFT10 (C3 grass) ! (18) ! Mixed Forest/woodland => PFT6 (Temperate Broadleaf Summergreen) ! (19) ! Interrupted Forest => PFT6 (Temperate Broadleaf Summergreen) ! (20) ! Water and Land Mixtures => none ! INTEGER(KIND=JPIM), INTENT(IN) :: KLON INTEGER(KIND=JPIM), INTENT(IN) :: KIDIA INTEGER(KIND=JPIM), INTENT(IN) :: KFDIA LOGICAL , INTENT(IN) :: LDLAND(:) INTEGER(KIND=JPIM), INTENT(IN) :: KVTYPE(:) TYPE(TAGF) , INTENT(IN) :: YDAGF INTEGER(KIND=JPIM), INTENT(IN) :: KSTEP REAL(KIND=JPRB) , INTENT(IN) :: dt_sechiba REAL(KIND=JPRB), INTENT(IN) :: swdown(:) ! Downwelling short wave flux ! @tex ($W m^{-2}$) @endtex REAL(KIND=JPRB), INTENT(IN) :: pb(:) ! Lowest level pressure (hPa) REAL(KIND=JPRB), INTENT(IN) :: qsurf(:) ! Near surface specific humidity ! @tex ($kg kg^{-1}$) @endtex REAL(KIND=JPRB), INTENT(IN) :: qair(:) ! Specific humidity at first atmospheric model layer ! @tex ($kg kg^{-1}$) @endtex REAL(KIND=JPRB), INTENT(IN) :: temp_air(:) ! Air temperature at first atmospheric model layer (K) REAL(KIND=JPRB), INTENT(IN) :: skin_temp(:) ! Surface temperature (K) !REAL(KIND=JPRB), INTENT(IN) :: rau(:) ! air density @tex ($kg m^{-3}$) @endtex !REAL(KIND=JPRB), INTENT(IN) :: u(:) ! Lowest level wind speed ! @tex ($m s^{-1}$) @endtex !REAL(KIND=JPRB), INTENT(IN) :: v(:) ! Lowest level wind speed ! @tex ($m s^{-1}$) @endtex !REAL(KIND=JPRB), INTENT(IN) :: q_cdrag(:) ! Surface drag coefficient (-) !REAL(KIND=JPRB), INTENT(IN) :: Ca(:) ! CO2 concentration inside the canopy REAL(KIND=JPRB), INTENT(IN) :: Ca ! CO2 concentration inside the canopy ! @tex ($\mu mol mol^{-1}$) @endtex REAL(KIND=JPRB), INTENT(IN) :: humrel(:) ! Soil moisture stress (0-1,unitless) !REAL(KIND=JPRB), INTENT(IN) :: veget(:) ! Coverage fraction of vegetation for each PFT ! depending on LAI (0-1, unitless) !REAL(KIND=JPRB), INTENT(IN) :: veget_max(:) ! Maximum vegetation fraction of each PFT inside ! the grid box (0-1, unitless) REAL(KIND=JPRB), INTENT(IN) :: lai(:) ! Leaf area index @tex ($m^2 m^{-2}$) @endtex ! @tex ($m s^{-1}$) @endtex !REAL(KIND=JPRB), INTENT(IN) :: qsintveg(:) ! Water on vegetation due to interception ! @tex ($kg m^{-2}$) @endte !REAL(KIND=JPRB), INTENT(IN) :: qsintmax(:) ! Maximum water on vegetation ! @tex ($kg m^{-2}$) @endtex !REAL(KIND=JPRB), INTENT(OUT) :: vbeta3(:,:) ! Beta for Transpiration (unitless) !REAL(KIND=JPRB), INTENT(OUT) :: vbeta3pot(:,:) ! Beta for Potential Transpiration !REAL(KIND=JPRB), INTENT(OUT) :: cimean(:,:) ! mean intercellular CO2 concentration ! @tex ($\mu mol mol^{-1}$) @endtex REAL(KIND=JPRB), INTENT(OUT) :: gsmean(:) ! mean stomatal conductance to CO2 (umol m-2 s-1) REAL(KIND=JPRB), INTENT(OUT) :: gpp(:) ! Assimilation (gC m^{-2} s^{-1}) !REAL(KIND=JPRB), INTENT(OUT) :: cim(:,:) ! Intercellular CO2 over nlai !* 0. LOCAL VARIABLES. ! ----- ---------- REAL(KIND=JPRB) :: rveget(KLON) ! stomatal resistance of vegetation ! @tex ($s m^{-1}$) @endtex REAL(KIND=JPRB) :: rstruct(KLON) ! structural resistance @tex ($s m^{-1}$) @endtex REAL(KIND=JPRB) :: vbeta23(KLON) ! Beta for fraction of wetted foliage that will ! transpire (unitless) INTEGER(KIND=JPIM) :: ICO2TYPE ! type of CO2 vegetation REAL(KIND=JPRB) :: vcmax(KLON) ! maximum rate of carboxylation ! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex INTEGER(KIND=JPIM) :: ji, jv, jl, limit_photo, JLEVEL ! indices (unitless) INTEGER(KIND=JPIM) :: info_limitphoto(KLON,21) REAL(KIND=JPRB) :: leaf_ci(KLON,20) ! intercellular CO2 concentration (ppm) ! ZCI = Leaf internal concentration of CO2 (kgCO2 kgAir-1) REAL(KIND=JPRB) :: leaf_ci_lowest(KLON) ! intercellular CO2 concentration at the lowest ! LAI level ! @tex ($\mu mol mol^{-1}$) @endtex INTEGER(KIND=JPIM) :: ilai(KLON) ! counter for loops on LAI levels (unitless) !REAL(KIND=JPRB) :: zqsvegrap(KLON) ! relative water quantity in the water ! interception reservoir (0-1,unitless) REAL(KIND=JPRB) :: speed ! wind speed @tex ($m s^{-1}$) @endtex ! Assimilation LOGICAL :: assimilate(KLON) ! where assimilation is to be calculated (unitless) LOGICAL :: calculate(KLON) ! where assimilation is to be calculated for in the PFTs loop (unitless) INTEGER(KIND=JPIM) :: nic,inic,icinic ! counter/indices (unitless) INTEGER(KIND=JPIM) :: index_calc(KLON) ! index (unitless) INTEGER(KIND=JPIM) :: nia,inia,nina,inina,iainia ! counter/indices (unitless) INTEGER(KIND=JPIM) :: index_assi(KLON),index_non_assi(KLON) ! indices (unitless) REAL(KIND=JPRB) :: vc2(KLON,21) ! rate of carboxylation (at a specific LAI level) ! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex REAL(KIND=JPRB) :: vj2(KLON,21) ! rate of Rubisco regeneration (at a specific LAI ! level) @tex ($\mu mol e- m^{-2} s^{-1}$) @endtex REAL(KIND=JPRB) :: assimi(KLON,21) ! assimilation (at a specific LAI level) ! @tex ($\mu mol m^{-2} s^{-1}$) @endtex ! (temporary variables) REAL(KIND=JPRB) :: gstop(KLON) ! stomatal conductance to H2O at topmost level ! @tex ($m s^{-1}$) @endtex REAL(KIND=JPRB) :: gs(KLON,21) ! stomatal conductance to CO2 ! @tex ($\mol m^{-2} s^{-1}$) @endtex REAL(KIND=JPRB) :: templeafci(KLON,20) REAL(KIND=JPRB) :: gamma_star(KLON) ! CO2 compensation point (ppm) ! @tex ($\mu mol mol^{-1}$) @endtex !REAL(KIND=JPRB) :: air_relhum(KLON) ! air relative humidity at 2m ! @tex ($kg kg^{-1}$) @endtex REAL(KIND=JPRB) :: VPD(KLON) ! Vapor Pressure Deficit (kPa) REAL(KIND=JPRB) :: water_lim(KLON) ! water limitation factor (0-1,unitless) REAL(KIND=JPRB) :: gstot(KLON) ! total stomatal conductance to H2O ! Final unit is ! @tex ($m s^{-1}$) @endtex REAL(KIND=JPRB) :: assimtot(KLON) ! total assimilation ! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex REAL(KIND=JPRB) :: Rdtot(KLON) ! Total Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 m−2 s−1) REAL(KIND=JPRB) :: leaf_gs_top(KLON) ! leaf stomatal conductance to H2O at topmost level ! @tex ($\mol H2O m^{-2} s^{-1}$) @endtex REAL(KIND=JPRB) :: laitab(21) ! tabulated LAI steps @tex ($m^2 m^{-2}$) @endtex REAL(KIND=JPRB) :: qsatt(KLON) ! surface saturated humidity at 2m (??) ! @tex ($g g^{-1}$) @endtex REAL(KIND=JPRB) :: light(20) ! fraction of light that gets through upper LAI ! levels (0-1,unitless) REAL(KIND=JPRB) :: T_Vcmax(KLON) ! Temperature dependance of Vcmax (unitless) REAL(KIND=JPRB) :: S_Vcmax_acclim_temp(KLON) ! Entropy term for Vcmax ! accounting for acclimation to temperature (J K-1 mol-1) REAL(KIND=JPRB) :: T_Jmax(KLON) ! Temperature dependance of Jmax REAL(KIND=JPRB) :: S_Jmax_acclim_temp(KLON) ! Entropy term for Jmax ! accounting for acclimation toxs temperature (J K-1 mol-1) REAL(KIND=JPRB) :: T_gm(KLON) ! Temperature dependance of gmw REAL(KIND=JPRB) :: T_Rd(KLON) ! Temperature dependance of Rd (unitless) REAL(KIND=JPRB) :: T_Kmc(KLON) ! Temperature dependance of KmC (unitless) REAL(KIND=JPRB) :: T_KmO(KLON) ! Temperature dependance of KmO (unitless) REAL(KIND=JPRB) :: T_Sco(KLON) ! Temperature dependance of Sco REAL(KIND=JPRB) :: T_gamma_star(KLON) ! Temperature dependance of gamma_star (unitless) REAL(KIND=JPRB) :: vc(KLON) ! Maximum rate of Rubisco activity-limited carboxylation (mumol CO2 m−2 s−1) REAL(KIND=JPRB) :: vj(KLON) ! Maximum rate of e- transport under saturated light (mumol CO2 m−2 s−1) REAL(KIND=JPRB) :: gm(KLON) ! Mesophyll diffusion conductance (molCO2 m−2 s−1 bar−1) REAL(KIND=JPRB) :: g0var(KLON) REAL(KIND=JPRB) :: Rd(KLON,21) ! Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 m−2 s−1) REAL(KIND=JPRB) :: Kmc(KLON) ! Michaelis–Menten constant of Rubisco for CO2 (mubar) REAL(KIND=JPRB) :: KmO(KLON) ! Michaelis–Menten constant of Rubisco for O2 (mubar) REAL(KIND=JPRB) :: Sco(KLON) ! Relative CO2 /O2 specificity factor for Rubisco (bar bar-1) REAL(KIND=JPRB) :: gb_co2(KLON) ! Boundary-layer conductance (molCO2 m−2 s−1 bar−1) REAL(KIND=JPRB) :: gb_h2o(KLON) ! Boundary-layer conductance (molH2O m−2 s−1 bar−1) REAL(KIND=JPRB) :: fvpd(KLON) ! Factor for describing the effect of leaf-to-air vapour difference on gs (-) REAL(KIND=JPRB) :: low_gamma_star(KLON) ! Half of the reciprocal of Sc/o (bar bar-1) REAL(KIND=JPRB) :: N_Vcmax ! Nitrogen level dependance of Vcmacx and Jmax REAL(KIND=JPRB) :: fcyc ! Fraction of electrons at PSI that follow cyclic transport around PSI (-) REAL(KIND=JPRB) :: z ! A lumped parameter (see Yin et al. 2009) ( mol mol-1) REAL(KIND=JPRB) :: Rm ! Day respiration in the mesophyll (umol CO2 m−2 s−1) REAL(KIND=JPRB) :: Cs_star ! Cs -based CO2 compensation point in the absence of Rd (ubar) REAL(KIND=JPRB) :: Iabs(KLON) ! Photon flux density absorbed by leaf photosynthetic pigments (umol photon m−2 s−1) REAL(KIND=JPRB) :: Jmax(KLON) ! Maximum value of J under saturated light (umol e− m−2 s−1) REAL(KIND=JPRB) :: JJ(KLON,21) ! Rate of e− transport (umol e− m−2 s−1) REAL(KIND=JPRB) :: J2 ! Rate of all e− transport through PSII (umol e− m−2 s−1) REAL(KIND=JPRB) :: VpJ2 ! e− transport-limited PEP carboxylation rate (umol CO2 m−2 s−1) REAL(KIND=JPRB) :: A_1, A_3 ! Lowest First and third roots of the analytical solution for a general cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 m−2 s−1) REAL(KIND=JPRB) :: A_1_tmp, A_3_tmp ! Temporary First and third roots of the analytical solution for a general cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 m−2 s−1) REAL(KIND=JPRB) :: Obs ! Bundle-sheath oxygen partial pressure (ubar) REAL(KIND=JPRB) :: Cc(KLON,21) ! Chloroplast CO2 partial pressure (ubar) REAL(KIND=JPRB) :: ci_star ! Ci-based CO2 compensation point in the absence of Rd (ubar) REAL(KIND=JPRB) :: a,b,c,d,m,f,j,g,h,i,l,p,q,r ! Variables used for solving the cubic equation (see Yin et al. (2009)) REAL(KIND=JPRB) :: QQ,UU,PSI,x1,x2,x3 ! Variables used for solving the cubic equation (see Yin et al. (2009)) REAL(KIND=JPRB) :: cresist ! coefficient for resistances (??) INTEGER(KIND=JPIM) :: KVTYPE_TEMPO(KLON) ! to deal with KVTYPE=0 !REAL(KIND=JPRB) :: laisum(KLON) ! when calculating cim over nlai !REAL(KIND=JPRB) :: wind(KLON) ! Wind (m s^{-1}) LOGICAL :: downregulation_co2 REAL(KIND=JPRB) :: leaf_temp(KLON) ! Leaf_temperature (?K) LOGICAL :: skin_temp_OK ! Logical to select leaf temperature between skin and surface air temperatures REAL(KIND=JPRB) :: ZHOOK_HANDLE IF (LHOOK) CALL DR_HOOK('FARQUHAR_MOD:FARQUHAR',0,ZHOOK_HANDLE) ASSOCIATE(max_temp=>YDAGF%max_temp, qsfrict=>YDAGF%qsfrict, & & nlai=>YDAGF%nlai, NPATH=>YDAGF%NPATH, zero=>YDAGF%zero, one=>YDAGF%one, two=>YDAGF%two, three=>YDAGF%three, & & four=>YDAGF%four, RR=>YDAGF%RR, min_wind=>YDAGF%min_wind, msmlr_h2o=>YDAGF%msmlr_h2o, & & downregulation_co2_baselevel=>YDAGF%downregulation_co2_baselevel, downregulation_co2_minimum=>YDAGF%downregulation_co2_minimum, & & pi=>YDAGF%pi, one_point_five=>YDAGF%one_point_five, nine=>YDAGF%nine, twenty_seven=>YDAGF%twenty_seven, & & fifty_four=>YDAGF%fifty_four, undef_sechiba=>YDAGF%undef_sechiba, undef=>YDAGF%undef, one_month=>YDAGF%one_month, a1=>YDAGF%a1, & & b1=>YDAGF%b1, arJV=>YDAGF%arJV, brJV=>YDAGF%brJV, aSJ=>YDAGF%aSJ, bSJ=>YDAGF%bSJ, aSV=>YDAGF%aSV, bSV=>YDAGF%bSV, & & D_Vcmax=>YDAGF%D_Vcmax, D_Jmax=>YDAGF%D_Jmax, E_gamma_star=>YDAGF%E_gamma_star, E_gm=>YDAGF%E_gm, S_gm=>YDAGF%S_gm, & & D_gm=>YDAGF%D_gm, E_Vcmax=>YDAGF%E_Vcmax, E_Jmax=>YDAGF%E_Jmax, E_KmC=>YDAGF%E_KmC, E_KmO=>YDAGF%E_KmO, E_Rd=>YDAGF%E_Rd, & & E_Sco=>YDAGF%E_Sco, g0=>YDAGF%g0, gamma_star25=>YDAGF%gamma_star25, gm25=>YDAGF%gm25, KmC25=>YDAGF%KmC25, KmO25=>YDAGF%KmO25, & & Sco25=>YDAGF%Sco25, stress_gm=>YDAGF%stress_gm, stress_gs=>YDAGF%stress_gs, stress_vcmax=>YDAGF%stress_vcmax, & & tphoto_max=>YDAGF%tphoto_max, tphoto_min=>YDAGF%tphoto_min, Vcmax25=>YDAGF%Vcmax25, & & downregulation_co2_coeff=>YDAGF%downregulation_co2_coeff, rstruct_const=>YDAGF%rstruct_const, CO2TYPE=>YDAGF%CO2TYPE, & & lai_level_depth=>YDAGF%lai_level_depth, laimax=>YDAGF%laimax, ext_coeff=>YDAGF%ext_coeff, gb_ref=>YDAGF%gb_ref, & & tp_00=>YDAGF%tp_00, pb_std=>YDAGF%pb_std, RG_to_PAR=>YDAGF%RG_to_PAR, W_to_mol=>YDAGF%W_to_mol, alpha_LL=>YDAGF%alpha_LL, & & theta=>YDAGF%theta, fpsir=>YDAGF%fpsir, fQ=>YDAGF%fQ, fpseudo=>YDAGF%fpseudo, h_protons=>YDAGF%h_protons, gbs=>YDAGF%gbs, & & Oi=>YDAGF%Oi, alpha=>YDAGF%alpha, kp=>YDAGF%kp, Tetens_1=>YDAGF%Tetens_1, Tetens_2=>YDAGF%Tetens_2, ref_temp=>YDAGF%ref_temp, & & mol_to_m_1=>YDAGF%mol_to_m_1, ratio_H2O_to_CO2=>YDAGF%ratio_H2O_to_CO2, N_vert_att=>YDAGF%N_vert_att, rveg_pft=>YDAGF%rveg_pft, & & tau_temp_air_month=>YDAGF%tau_temp_air_month) ! ! initialization of local variables ! downregulation_co2 = .TRUE. ! Set to .TRUE. if you want CO2 downregulation. skin_temp_OK = .TRUE. ! Leaf temperature is equal to skin temperature (else to surface air temperature) IF (FIRST_CALL) THEN !CALL qsfrict_init !TO BE RESTARTED ALLOCATE (temp_growth(KLON)) temp_growth(:) = temp_air(:) - tp_00 FIRST_CALL=.FALSE. ENDIF !temp_growth is updated only once per time step, not each time farquhar is called for the differene vegetation types. ! We use temp_air (Kattge & Knorr, 2007: "the plants’ average ambient temperature"). IF (KSTEP .NE. KSTEP_SAV) THEN IF (KSTEP .LT. one_month/dt_sechiba) THEN !Simple mean while lower than one month temp_growth(:) = ( temp_growth(:) * KSTEP + temp_air(:) - tp_00) / (KSTEP+1) ELSE !Weighted average temp_growth(:) = ( temp_growth(:) * ( tau_temp_air_month - dt_sechiba ) + (temp_air(:) - tp_00) * dt_sechiba ) / tau_temp_air_month ENDIF WHERE ( ABS(temp_growth(:) + tp_00) .LT. EPSILON(zero) ) temp_growth(:) = -tp_00 ENDWHERE ENDIF IF (skin_temp_OK) THEN leaf_temp(:)=skin_temp(:) ELSE leaf_temp(:)=temp_air(:) ENDIF !wind(:) = SQRT (u(:)*u(:) + v(:)*v(:)) !! 1. Calculate the different coefficients !IF (.NOT.ldq_cdrag_from_gcm) THEN ! Case 1a) ! CALL diffuco_aero (KLON, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, & ! qsurf, qair, q_cdrag) !ENDIF ! Case 1b) !CALL diffuco_raerod (KLON, u, v, q_cdrag, raero) !! 2. Make an estimation of the saturated humidity at the surface !CALL qsatcalc (KLON, temp_sol, pb, qsatt) !! 4. Calculate the beta coefficient for interception !CALL diffuco_inter (KLON, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, & ! & qsintveg, qsintmax, rstruct, vbeta2, vbeta23) ! 1. Preliminary calculations !_ ================================================================================================================================ !cim(:,:)=zero ! ! 1.1 Calculate LAI steps ! The integration at the canopy level is done over nlai fixed levels. DO JLEVEL = 1, nlai+1 laitab(JLEVEL) = laimax*(EXP(lai_level_depth*REAL(JLEVEL-1,JPRB))-1.)/(EXP(lai_level_depth*REAL(nlai,JPRB))-one) ENDDO ! ! 1.2 Calculate light fraction for each LAI step ! The available light follows a simple Beer extinction law. ! The extinction coefficients (ext_coef) are PFT-dependant constants. DO JLEVEL = 1, nlai light(JLEVEL) = exp( -ext_coeff*laitab(JLEVEL) ) ENDDO !Here we set a protection for when KVTYPE(JL) is zero. We set the KVTYPE_TEMPO array to one in this case to be able to index the arrays that depend on KVTPE !and start at 1. A protection against KVTYPE being zero is explicitely done is the loops below to ensure that there are no caculations done in this case. WHERE(KVTYPE(:) .GT. 0) KVTYPE_TEMPO(:)=KVTYPE(:) ELSEWHERE KVTYPE_TEMPO(:)=1 ENDWHERE IF (downregulation_co2) THEN vcmax(:) = Vcmax25(KVTYPE_TEMPO(:))*(one-downregulation_co2_coeff(KVTYPE_TEMPO(:)) * & (MAX(Ca,downregulation_co2_minimum)-downregulation_co2_baselevel) / & (MAX(Ca,downregulation_co2_minimum)+20._JPRB)) ELSE vcmax(:) = Vcmax25(KVTYPE_TEMPO(:)) ENDIF ! ! 1.3 Estimate relative humidity of air (for calculation of the stomatal conductance). ! CALL qsatcalc (KIDIA, KFDIA, KLON, YDAGF, leaf_temp, pb, qsatt) !air_relhum(:) = & ! ( qair(:) * pb(:) / (Tetens_1+qair(:)* Tetens_2) ) / & ! ( qsatt(:)*pb(:) / (Tetens_1+qsatt(:)*Tetens_2 ) ) VPD(:) = ( qsatt(:)*pb(:) / (Tetens_1+qsatt(:)*Tetens_2 ) ) & - ( qair(:) * pb(:) / (Tetens_1+qair(:)* Tetens_2) ) ! VPD is needed in kPa VPD(:) = VPD(:)/10. ! ! 2. beta coefficient for vegetation transpiration ! !TO BE DEALT WITH FOR BARE SOIL !rstruct(:,1) = rstruct_const(1) rveget(:) = undef_sechiba ! !vbeta3(:,:) = zero !vbeta3pot(:,:) = zero gsmean(:) = zero gpp(:) = zero ! ! cimean(:,1) = Ca(:) ! ! 2. Loop over vegetation types ! !DO JVT=0,KVTYPES gamma_star(:)=zero Kmo(:)=zero Kmc(:)=zero gm(:)=zero g0var(:) =zero Cc(:,:)=zero Vc2(:,:)=zero JJ(:,:)=zero info_limitphoto(:,:)=0 gs(:,:)=zero templeafci(:,:)=zero assimi(:,:)=zero Rd(:,:)=zero ! ! 2.1 Initializations ! ! ! beta coefficient for vegetation transpiration ! rstruct(:) = rstruct_const(KVTYPE_TEMPO(:)) !cimean(:,JVT) = Ca(:) ! !! mask that contains points where there is photosynthesis !! For the sake of vectorisation [DISPENSABLE], computations are done only for convenient points. !! nia is the number of points where the assimilation is calculated and nina the number of points where photosynthesis is not !! calculated (based on criteria on minimum or maximum values on LAI, vegetation fraction, shortwave incoming radiation, !! temperature and relative humidity). !! For the points where assimilation is not calculated, variables are initialized to specific values. !! The assimilate(KLON) array contains the logical value (TRUE/FALSE) relative to this photosynthesis calculation. !! The index_assi(KLON) array indexes the nia points with assimilation, whereas the index_no_assi(KLON) array indexes !! the nina points with no assimilation. nia=0 nina=0 ! DO JL=KIDIA,KFDIA IF ( (LDLAND(JL)) .AND. ( lai(JL) .GT. 0.01 ) .AND. (KVTYPE(JL) .GT. 0)) THEN IF ( ( swdown(JL) .GT. EPSILON(1._JPRM) ) .AND. & ( humrel(JL) .GT. EPSILON(1._JPRM) ) .AND. & ( temp_growth(JL) .GT. tphoto_min ) .AND. & ( temp_growth(JL) .LT. tphoto_max ) ) THEN ! assimilate(JL) = .TRUE. nia=nia+1 index_assi(nia)=JL ! ELSE ! assimilate(JL) = .FALSE. nina=nina+1 index_non_assi(nina)=JL ! ENDIF ELSE ! assimilate(JL) = .FALSE. nina=nina+1 index_non_assi(nina)=JL ! ENDIF ! ENDDO ! gstot(:) = zero gstop(:) = zero assimtot(:) = zero Rdtot(:)=zero leaf_gs_top(:) = zero ! !zqsvegrap(:) = zero !WHERE (qsintmax(:) .GT. min_sechiba) !! relative water quantity in the water interception reservoir ! zqsvegrap(:) = MAX(zero, qsintveg(:) / qsintmax(:)) !ENDWHERE ! !! Calculates the water limitation factor. water_lim(:) = humrel(:) ! give a default value of ci for all pixel that do not assimilate DO JLEVEL=1,nlai DO inina=1,nina leaf_ci(index_non_assi(inina),JLEVEL) = Ca ENDDO ENDDO ! ilai(:) = 1 ! ! Here is the calculation of assimilation and stomatal conductance ! based on the work of Farquahr, von Caemmerer and Berry (FvCB model) ! as described in Yin et al. 2009 ! Yin et al. developed a extended version of the FvCB model for C4 plants ! and proposed an analytical solution for both photosynthesis pathways (C3 and C4) ! Photosynthetic parameters used are those reported in Yin et al. ! Except For Vcmax25, relationships between Vcmax25 and Jmax25 for which we use ! Medlyn et al. (2002) and Kattge & Knorr (2007) ! Because these 2 references do not consider mesophyll conductance, we neglect this term ! in the formulations developed by Yin et al. ! Consequently, gm (the mesophyll conductance) tends to the infinite ! This is of importance because as stated by Kattge & Knorr and Medlyn et al., ! values of Vcmax and Jmax derived with different model parametrizations are not ! directly comparable and the published values of Vcmax and Jmax had to be standardized ! to one consistent formulation and parametrization ! See eq. 6 of Yin et al. (2009) ! Parametrization of Medlyn et al. (2002) - from Bernacchi et al. (2001) T_KmC(:) = Arrhenius(KLON,leaf_temp,ref_temp,E_KmC,YDAGF) T_KmO(:) = Arrhenius(KLON,leaf_temp,ref_temp,E_KmO,YDAGF) T_Sco(:) = Arrhenius(KLON,leaf_temp,ref_temp,E_Sco,YDAGF) T_gamma_star(:) = Arrhenius(KLON,leaf_temp,ref_temp,E_gamma_star,YDAGF) ! Parametrization of Yin et al. (2009) - from Bernacchi et al. (2001) T_Rd(:) = Arrhenius(KLON,leaf_temp,ref_temp,E_Rd,YDAGF) ! For C3 plants, we assume that the Entropy term for Vcmax and Jmax ! acclimates to temperature as shown by Kattge & Knorr (2007) - Eq. 9 and 10 ! and that Jmax and Vcmax respond to temperature following a modified Arrhenius function ! (with a decrease of these parameters for high temperature) as in Medlyn et al. (2002) ! and Kattge & Knorr (2007). ! In Yin et al. (2009), temperature dependance to Vcmax is based only on a Arrhenius function ! Concerning this apparent unconsistency, have a look to the section 'Limitation of ! Photosynthesis by gm' of Bernacchi (2002) that may provide an explanation ! Growth temperature tested by Kattge & Knorr range from 11 to 35°C ! So, we limit the relationship between these lower and upper limits WHERE(CO2TYPE(KVTYPE_TEMPO(:)) .EQ. one) !C3 S_Jmax_acclim_temp(:) = aSJ(1) + bSJ(1) * MAX(11., MIN(temp_growth(:),35.)) T_Jmax(:) = Arrhenius_modified(KLON,leaf_temp,ref_temp,E_Jmax(1),D_Jmax(1),S_Jmax_acclim_temp,YDAGF) S_Vcmax_acclim_temp(:) = aSV(1) + bSV(1) * MAX(11., MIN(temp_growth(:),35.)) T_Vcmax(:) = Arrhenius_modified(KLON,leaf_temp,ref_temp,E_Vcmax(1),D_Vcmax(1),S_Vcmax_acclim_temp,YDAGF) ELSEWHERE(CO2TYPE(KVTYPE_TEMPO(:)) .EQ. two) !C4 S_Jmax_acclim_temp(:) = aSJ(2) + bSJ(2) * MAX(11., MIN(temp_growth(:),35.)) T_Jmax(:) = Arrhenius_modified(KLON,leaf_temp,ref_temp,E_Jmax(2),D_Jmax(2),S_Jmax_acclim_temp,YDAGF) S_Vcmax_acclim_temp(:) = aSV(2) + bSV(2) * MAX(11., MIN(temp_growth(:),35.)) T_Vcmax(:) = Arrhenius_modified(KLON,leaf_temp,ref_temp,E_Vcmax(2),D_Vcmax(2),S_Vcmax_acclim_temp,YDAGF) ENDWHERE vc(:) = vcmax(:) * T_Vcmax(:) ! As shown by Kattge & Knorr (2007), we make use ! of Jmax25/Vcmax25 ratio (rJV) that acclimates to temperature for C3 plants ! rJV is written as a function of the growth temperature ! rJV = arJV + brJV * T_month ! See eq. 10 of Kattge & Knorr (2007) ! and Table 3 for Values of arJV anf brJV ! Growth temperature is monthly temperature (expressed in °C) - See first paragraph of ! section Methods/Data of Kattge & Knorr WHERE(CO2TYPE(KVTYPE_TEMPO(:)) .EQ. one) !C3 vj(:) = ( arJV(1) + brJV(1) * MAX(11., MIN(temp_growth(:),35.)) ) * vcmax(:) * T_Jmax(:) T_gm(:) = Arrhenius_modified(KLON,leaf_temp,ref_temp,E_gm(1),D_gm(1),S_gm(1),YDAGF) gm(:) = gm25(1) * T_gm(:) * MAX(1-stress_gm, water_lim(:)) g0var(:) = g0(1)* MAX(1-stress_gs, water_lim(:)) KmC(:)=KmC25(1)*T_KmC(:) KmO(:)=KmO25(1)*T_KmO(:) Sco(:)=Sco25(1)*T_sco(:) ELSEWHERE(CO2TYPE(KVTYPE_TEMPO(:)) .EQ. two) !C4 vj(:) = ( arJV(2) + brJV(2) * MAX(11., MIN(temp_growth(:),35.)) ) * vcmax(:) * T_Jmax(:) T_gm(:) = Arrhenius_modified(KLON,leaf_temp,ref_temp,E_gm(2),D_gm(2),S_gm(2),YDAGF) !gm not neded for C4 plants !gm(:) = gm25(2) * T_gm(:) * MAX(1-stress_gm, water_lim(:)) g0var(:) = g0(2)* MAX(1-stress_gs, water_lim(:)) KmC(:)=KmC25(2)*T_KmC(:) KmO(:)=KmO25(2)*T_KmO(:) Sco(:)=Sco25(2)*T_sco(:) ELSEWHERE Sco(:)=Sco25(2) ENDWHERE gamma_star(:) = gamma_star25*T_gamma_star(:) ! low_gamma_star is defined by Yin et al. (2009) ! as the half of the reciprocal of Sco - See Table 2 low_gamma_star(:) = 0.5 / Sco(:) ! VPD expressed in kPa ! Note : MIN(1.-min_sechiba,MAX(min_sechiba,(a1(JVT) - b1(JVT) * VPD(:)))) is always between 0-1 not including 0 and 1 WHERE(CO2TYPE(KVTYPE_TEMPO(:)) .EQ. one) !C3 fvpd(:) = 1. / ( 1. / MIN(1.-EPSILON(1._JPRM),MAX(EPSILON(1._JPRM),(a1(1) - b1(1) * VPD(:)))) - 1. ) & * MAX(1-stress_gs, water_lim(:)) ELSEWHERE(CO2TYPE(KVTYPE_TEMPO(:)) .EQ. two) !C4 fvpd(:) = 1. / ( 1. / MIN(1.-EPSILON(1._JPRM),MAX(EPSILON(1._JPRM),(a1(2) - b1(2) * VPD(:)))) - 1. ) & * MAX(1-stress_gs, water_lim(:)) ENDWHERE ! leaf boundary layer conductance ! conversion from a conductance in (m s-1) to (mol H2O m-2 s-1) ! from Pearcy et al. (1991, see below) gb_h2o(:) = gb_ref * 44.6 * (tp_00/temp_air(:)) * (pb(:)/pb_std) ! conversion from (mol H2O m-2 s-1) to (mol CO2 m-2 s-1) gb_co2(:) = gb_h2o(:) / ratio_H2O_to_CO2 ! ! 2.4 Loop over LAI discretized levels to estimate assimilation and conductance ! !! The calculate(KLON) array is of type logical to indicate wether we have to sum over this LAI fixed level (the LAI of !! 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 !! corresponding point is indexed in the index_calc array. JJ(:,:)=zero vc2(:,:)=zero vj2(:,:)=zero Cc(:,:)=zero gs(:,:)=zero assimi(:,:)=zero Rd(:,:)=zero DO JLEVEL = 1, nlai ! nic=0 calculate(:) = .FALSE. ! IF (nia .GT. 0) then DO inia=1,nia calculate(index_assi(inia)) = (laitab(JLEVEL) .LE. lai(index_assi(inia)) ) IF ( calculate(index_assi(inia)) ) THEN nic=nic+1 index_calc(nic)=index_assi(inia) ENDIF ENDDO ENDIF ! ! ! 2.4.1 Vmax is scaled into the canopy due to reduction of nitrogen !! (Johnson and Thornley,1984). ! N_Vcmax = ( one - N_vert_att * ( one - light(JLEVEL) ) ) ! vc2(:,JLEVEL) = vc(:) * N_Vcmax * MAX(1-stress_vcmax, water_lim(:)) vj2(:,JLEVEL) = vj(:) * N_Vcmax * MAX(1-stress_vcmax, water_lim(:)) ! see Comment in legend of Fig. 6 of Yin et al. (2009) ! Rd25 is assumed to equal 0.01 Vcmax25 Rd(:,JLEVEL) = vcmax(:) * N_Vcmax * 0.01 * T_Rd(:) * MAX(1-stress_vcmax, water_lim(:)) Iabs(:)=swdown(:)*W_to_mol*RG_to_PAR*ext_coeff*light(JLEVEL) ! eq. 4 of Yin et al (2009) Jmax(:)=vj2(:,JLEVEL) JJ(:,JLEVEL) = ( alpha_LL * Iabs(:) + Jmax(:) - sqrt((alpha_LL * Iabs(:) + Jmax(:) )**2. & - 4 * theta * Jmax(:) * alpha_LL * Iabs(:)) ) & / ( 2 * theta) ! !IF ( CO2TYPE(KVTYPE)) .EQ. one ) THEN ! ! 2.4.2 Assimilation for C4 plants (Collatz et al., 1992) ! IF (nic .GT. 0) THEN DO inic=1,nic ! Analytical resolution of the Assimilation based Yin et al. (2009) icinic=index_calc(inic) IF (CO2TYPE(KVTYPE(icinic)) .EQ. two) THEN !C4 ! Eq. 28 of Yin et al. (2009) fcyc= 1. - ( 4.*(1.-fpsir)*(1.+fQ) + 3.*h_protons*fpseudo ) / & ( 3.*h_protons - 4.*(1.-fpsir)) ! See paragraph after eq. (20b) of Yin et al. Rm=Rd(icinic,JLEVEL)/2. ! We assume that cs_star equals ci_star (see Comment in legend of Fig. 6 of Yin et al. (2009) ! Equation 26 of Yin et al. (2009) Cs_star = (gbs * low_gamma_star(icinic) * Oi - & ( 1. + low_gamma_star(icinic) * alpha / 0.047) * Rd(icinic,JLEVEL) + Rm ) & / ( gbs + kp ) ! eq. 11 of Yin et al (2009) J2 = JJ(icinic,JLEVEL) / ( 1. - fpseudo / ( 1. - fcyc ) ) ! Equation right after eq. (20d) of Yin et al. (2009) z = ( 2. + fQ - fcyc ) / ( h_protons * (1. - fcyc )) VpJ2 = fpsir * J2 * z / 2. A_3=9999. ! See eq. right after eq. 18 of Yin et al. (2009) DO limit_photo=1,2 ! Is Vc limiting the Assimilation IF ( limit_photo .EQ. 1 ) THEN a = 1. + kp / gbs b = 0. x1 = Vc2(icinic,JLEVEL) x2 = KmC(icinic)/KmO(icinic) x3 = KmC(icinic) ! Is J limiting the Assimilation ELSE a = 1. b = VpJ2 x1 = (1.- fpsir) * J2 * z / 3. x2 = 7. * low_gamma_star(icinic) / 3. x3 = 0. ENDIF m=fvpd(icinic)-g0var(icinic)/gb_co2(icinic) d=g0var(icinic)*(Ca-Cs_star) + fvpd(icinic)*Rd(icinic,JLEVEL) f=(b-Rm-low_gamma_star(icinic)*Oi*gbs)*x1*d + a*gbs*x1*Ca*d j=(b-Rm+gbs*x3 + x2*gbs*Oi)*m + (alpha*x2/0.047-1.)*d & + a*gbs*(Ca*m - d/gb_co2(icinic) - (Ca - Cs_star )) g=(b-Rm-low_gamma_star(icinic)*Oi*gbs)*x1*m - (alpha*low_gamma_star(icinic)/0.047+1.)*x1*d & + a*gbs*x1*(Ca*m - d/gb_co2(icinic) - (Ca-Cs_star )) h=-((alpha*low_gamma_star(icinic)/0.047+1.)*x1*m + (a*gbs*x1*(m-1.))/gb_co2(icinic) ) i= ( b-Rm + gbs*x3 + x2*gbs*Oi )*d + a*gbs*Ca*d l= ( alpha*x2/0.047 - 1.)*m - (a*gbs*(m-1.))/gb_co2(icinic) p = (j-(h-l*Rd(icinic,JLEVEL))) / l q = (i+j*Rd(icinic,JLEVEL)-g) / l r = -(f-i*Rd(icinic,JLEVEL)) / l ! See Yin et al. (2009) and Baldocchi (1994) QQ = ( (p**two) - three * q) / nine UU = ( two* (p**three) - nine *p*q + twenty_seven *r) /fifty_four IF ( (QQ .GE. zero) .AND. (ABS(UU/(ABS(QQ)**one_point_five) ) .LE. one) ) THEN PSI = ACOS(MAX(-1._JPRB,MIN(UU/(QQ**one_point_five),one))) A_3_tmp = -two * SQRT(QQ) * COS(( PSI + four * pi)/three ) - p / three IF (( A_3_tmp .LT. A_3 )) THEN A_3 = A_3_tmp info_limitphoto(icinic,JLEVEL)=2 ELSE ! In case, J is not limiting the assimilation ! we have to re-initialise a, b, x1, x2 and x3 values ! in agreement with a Vc-limited assimilation a = one + kp / gbs b = zero x1 = Vc2(icinic,JLEVEL) x2 = KmC(icinic)/KmO(icinic) x3 = KmC(icinic) info_limitphoto(icinic,JLEVEL)=1 ENDIF ENDIF IF ( ( A_3 .EQ. 9999. ) .OR. ( A_3 .LT. (-Rd(icinic,JLEVEL)) ) ) THEN A_3 = -Rd(icinic,JLEVEL) ENDIF assimi(icinic,JLEVEL) = A_3 IF ( ABS( assimi(icinic,JLEVEL) + Rd(icinic,JLEVEL) ) .LT. EPSILON(1._JPRM) ) THEN gs(icinic,JLEVEL) = g0var(icinic) !leaf_ci keeps its initial value (Ca). ELSE ! Eq. 24 of Yin et al. (2009) Obs = ( alpha * assimi(icinic,JLEVEL) ) / ( 0.047 * gbs ) + Oi ! Eq. 23 of Yin et al. (2009) Cc(icinic,JLEVEL) = ( ( assimi(icinic,JLEVEL) + Rd(icinic,JLEVEL) ) * ( x2 * Obs + x3 ) + low_gamma_star(icinic) & * Obs * x1 ) & / MAX(EPSILON(1._JPRM), x1 - ( assimi(icinic,JLEVEL) + Rd(icinic,JLEVEL) )) ! Eq. 22 of Yin et al. (2009) leaf_ci(icinic,JLEVEL) = ( Cc(icinic,JLEVEL) - ( b - assimi(icinic,JLEVEL) - Rm ) / gbs ) / a ! Eq. 25 of Yin et al. (2009) ! It should be Cs instead of Ca but it seems that ! other equations in Appendix C make use of Ca gs(icinic,JLEVEL) = g0var(icinic) + ( assimi(icinic,JLEVEL) + Rd(icinic,JLEVEL) ) / & ( Ca - Cs_star ) * fvpd(icinic) ENDIF ENDDO !end limit_photo !ENDDO !end icinic !ENDIF !ELSE IF (CO2TYPE(KVEG(KVTYPE)) .EQ. two) THEN ELSE IF (CO2TYPE(KVTYPE(icinic)) .EQ. one) THEN ! ! 2.4.3 Assimilation for C3 plants (Farqhuar et al., 1980) ! !IF (nic .GT. 0) THEN !DO inic=1,nic !icinic=index_calc(inic) A_1=9999. ! See eq. right after eq. 18 of Yin et al. (2009) DO limit_photo=1,2 ! Is Vc limiting the Assimilation IF ( limit_photo .EQ. 1 ) THEN x1 = vc2(icinic,JLEVEL) ! It should be O not Oi (comment from Vuichard) x2 = KmC(icinic) * ( 1. + 2*gamma_star(icinic)*Sco(icinic) / KmO(icinic) ) ! Is J limiting the Assimilation ELSE x1 = JJ(icinic,JLEVEL)/4. x2 = two * gamma_star(icinic) ENDIF ! See Appendix B of Yin et al. (2009) a = g0var(icinic) * ( x2 + gamma_star(icinic) ) + & ( g0var(icinic) / gm(icinic) + fvpd(icinic) ) * ( x1 - Rd(icinic,JLEVEL) ) b = Ca * ( x1 - Rd(icinic,JLEVEL) ) - gamma_star(icinic) * x1 - Rd(icinic,JLEVEL) * x2 c = Ca + x2 + ( 1./gm(icinic) + 1./gb_co2(icinic) ) * ( x1 - Rd(icinic,JLEVEL) ) d = x2 + gamma_star(icinic) + ( x1 - Rd(icinic,JLEVEL) ) / gm(icinic) m = 1./gm(icinic) + ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * ( 1./gm(icinic) + 1./gb_co2(icinic) ) p = - ( d + (x1 - Rd(icinic,JLEVEL) ) / gm(icinic) + a * ( 1./gm(icinic) + 1./gb_co2(icinic) ) + & ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * c ) / m q = ( d * ( x1 - Rd(icinic,JLEVEL) ) + a*c + ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * b ) / m r = - a * b / m ! See Yin et al. (2009) QQ = ( (p**two) - three * q) / nine UU = ( two* (p**three) - nine *p*q + twenty_seven *r) /fifty_four IF ( (QQ .GE. zero) .AND. (ABS(UU/(QQ**one_point_five) ) .LE. one) ) THEN PSI = ACOS(UU/(QQ**one_point_five)) A_1_tmp = -two * SQRT(QQ) * COS( PSI / three ) - p / three IF (( A_1_tmp .LT. A_1 )) THEN A_1 = A_1_tmp info_limitphoto(icinic,JLEVEL)=2. ELSE ! In case, J is not limiting the assimilation ! we have to re-initialise x1 and x2 values ! in agreement with a Vc-limited assimilation x1 = vc2(icinic,JLEVEL) ! It should be O not Oi (comment from Vuichard) x2 = KmC(icinic) * ( one + 2*gamma_star(icinic)*Sco(icinic) / KmO(icinic) ) info_limitphoto(icinic,JLEVEL)=1 ENDIF ENDIF ENDDO !limi_photo IF ( (A_1 .EQ. 9999.) .OR. ( A_1 .LT. (-Rd(icinic,JLEVEL)) ) ) THEN A_1 = -Rd(icinic,JLEVEL) ENDIF assimi(icinic,JLEVEL) = A_1 IF ( ABS( assimi(icinic,JLEVEL) + Rd(icinic,JLEVEL) ) .LT. EPSILON(1._JPRM) ) THEN gs(icinic,JLEVEL) = g0var(icinic) ELSE ! Eq. 18 of Yin et al. (2009) Cc(icinic,JLEVEL) = ( gamma_star(icinic) * x1 + ( assimi(icinic,JLEVEL) + Rd(icinic,JLEVEL) ) * x2 ) & / MAX( EPSILON(1._JPRM), x1 - ( assimi(icinic,JLEVEL) + Rd(icinic,JLEVEL) ) ) ! Eq. 17 of Yin et al. (2009) leaf_ci(icinic,JLEVEL) = Cc(icinic,JLEVEL) + assimi(icinic,JLEVEL) / gm(icinic) ! See eq. right after eq. 15 of Yin et al. (2009) ci_star = gamma_star(icinic) - Rd(icinic,JLEVEL) / gm(icinic) ! ! Eq. 15 of Yin et al. (2009) gs(icinic,JLEVEL) = g0var(icinic) + ( assimi(icinic,JLEVEL) + Rd(icinic,JLEVEL) ) / ( leaf_ci(icinic,JLEVEL) & - ci_star ) * fvpd(icinic) ENDIF ELSE !STOP "Wrong ICO2TYPE in FARQUHAR" write (*,*) "Wrong ICO2TYPE in FARQUHAR" assimi(icinic,JLEVEL) = zero gs(icinic,JLEVEL) = g0var(icinic) ENDIF ENDDO !icinic ENDIF !nic ! IF (nic .GT. 0) THEN ! DO inic=1,nic ! !! 2.4.4 Estimatation of the stomatal conductance (Ball et al., 1987). ! icinic=index_calc(inic) ! ! keep stomatal conductance of topmost level ! IF ( JLEVEL .EQ. 1 ) THEN leaf_gs_top(icinic) = gs(icinic,JLEVEL) ! ENDIF ! !! 2.4.5 Integration at the canopy level ! total assimilation and conductance assimtot(icinic) = assimtot(icinic) + & assimi(icinic,JLEVEL) * (laitab(JLEVEL+1)-laitab(JLEVEL)) !Rdtot(icinic) = Rdtot(icinic) + & ! Rd(icinic,JLEVEL) * (laitab(JLEVEL+1)-laitab(JLEVEL)) gstot(icinic) = gstot(icinic) + & gs(icinic,JLEVEL) * (laitab(JLEVEL+1)-laitab(JLEVEL)) ! ilai(icinic) = JLEVEL ! ENDDO ! ENDIF ENDDO ! loop over LAI steps !! Calculated intercellular CO2 over nlai needed for the chemistry module !cim(:,JVT)=zero !laisum(:)=zero !DO JLEVEL=1,nlai ! WHERE (laitab(JLEVEL) .LE. lai(:,JVT) ) ! cim(:,JVT)= cim(:,JVT)+leaf_ci(:,JVT,JLEVEL)*(laitab(JLEVEL+1)-laitab(JLEVEL)) ! laisum(:)=laisum(:)+ (laitab(JLEVEL+1)-laitab(JLEVEL)) ! ENDWHERE !ENDDO !WHERE (laisum(:)>zero) ! cim(:,JVT)= cim(:,JVT)/laisum(:) !ENDWHERE ! !! 2.5 Calculate resistances ! IF (nia .GT. 0) THEN ! DO inia=1,nia ! iainia=index_assi(inia) !! Mean stomatal conductance for CO2 (mol m-2 s-1) gsmean(iainia) = gstot(iainia) ! ! cimean is the "mean ci" calculated in such a way that assimilation ! calculated in enerbil is equivalent to assimtot ! !IF ( ABS(gsmean(iainia,JVT)-g0var(iainia)*laisum(iainia)) .GT. min_sechiba) THEN ! cimean(iainia,JVT) = (fvpd(iainia)*(assimtot(iainia)+Rdtot(iainia))) /& ! (gsmean(iainia,JVT)-g0var(iainia)*laisum(iainia)) + gamma_star(iainia) !ELSE ! cimean(iainia,JVT) = gamma_star(iainia) !ENDIF ! conversion from umol m-2 (PFT) s-1 to kgCO2 m-2 (mesh area) tstep-1 !gpp(iainia,JVT) = assimtot(iainia)*12e-6*veget_max(iainia,JVT)*dt_sechiba !gpp(iainia) = assimtot(iainia)*12e-6 gpp(iainia) = assimtot(iainia)*44.e-9_JPRB !write(*,*) 'gpp(iainia)=',gpp(iainia) ! ! conversion from mol/m^2/s to m/s ! ! As in Pearcy, Schulze and Zimmermann ! Measurement of transpiration and leaf conductance ! Chapter 8 of Plant Physiological Ecology ! Field methods and instrumentation, 1991 ! Editors: ! ! Robert W. Pearcy, ! James R. Ehleringer, ! Harold A. Mooney, ! Philip W. Rundel ! ! ISBN: 978-0-412-40730-7 (Print) 978-94-010-9013-1 (Online) gstot(iainia) = mol_to_m_1 *(temp_air(iainia)/tp_00)*& (pb_std/pb(iainia))*gstot(iainia)*ratio_H2O_to_CO2 gstop(iainia) = mol_to_m_1 * (temp_air(iainia)/tp_00)*& (pb_std/pb(iainia))*leaf_gs_top(iainia)*ratio_H2O_to_CO2*& laitab(ilai(iainia)+1) ! IF (gstop(iainia) .LT. EPSILON(1._JPRM)) THEN rveget(iainia) = undef_sechiba ELSE rveget(iainia) = one/gstop(iainia) ENDIF ! ! ! rstruct is the difference between rtot (=1./gstot) and rveget ! ! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif IF (gstot(iainia) .LT. EPSILON(1._JPRM)) THEN rstruct(iainia) = undef_sechiba ELSE rstruct(iainia) = one/gstot(iainia) - & rveget(iainia) rstruct(iainia) = MAX( one/gstot(iainia) - & rveget(iainia), EPSILON(1._JPRM)) ENDIF ! ! !! wind is a global variable of the diffuco module. !speed = MAX(min_wind, wind(iainia)) ! ! beta for transpiration ! ! Corrections Nathalie - 28 March 2006 - on advices of Fred Hourdin !! Introduction of a potentiometer rveg_pft to settle the rveg+rstruct sum problem in the coupled mode. !! rveg_pft=1 in the offline mode. rveg_pft is a global variable declared in the diffuco module. !vbeta3(iainia,JVT) = veget_max(iainia,JVT) * & ! (one - zqsvegrap(iainia)) * & ! (one / (one + speed * q_cdrag(iainia) * (rveget(iainia,JVT) + & ! rstruct(iainia,JVT)))) !! Global resistance of the canopy to evaporation !cresist=(one / (one + speed * q_cdrag(iainia) * & ! veget(iainia,JVT)/veget_max(iainia,JVT) * & ! (rveg_pft*(rveget(iainia,JVT) + rstruct(iainia,JVT))))) !IF ( humrel(iainia,JVT) >= min_sechiba ) THEN ! vbeta3(iainia,JVT) = veget(iainia,JVT) * & ! (one - zqsvegrap(iainia)) * cresist + & ! MIN( vbeta23(iainia,JVT), veget(iainia,JVT) * & ! zqsvegrap(iainia) * cresist ) !ELSE ! Because of a minimum conductance g0, vbeta3 cannot be zero even if humrel=0 ! in the above equation. ! Here, we force transpiration to be zero when the soil cannot deliver it ! vbeta3(iainia,JVT) = zero !END IF ! vbeta3pot for computation of potential transpiration (needed for irrigation) !vbeta3pot(iainia,JVT) = MAX(zero, veget(iainia,JVT) * cresist) ! ! ENDDO ! ENDIF KSTEP_SAV=KSTEP ! !END DO ! loop over vegetation types ! END ASSOCIATE IF (LHOOK) CALL DR_HOOK('FARQUHAR_MOD:FARQUHAR',1,ZHOOK_HANDLE) END SUBROUTINE FARQUHAR !! ================================================================================================================================ !! SUBROUTINE : diffuco_aero !! !>\BRIEF This module first calculates the surface drag !! coefficient, for cases in which the surface drag coefficient is NOT provided by the coupled !! atmospheric model LMDZ or when the flag ldq_cdrag_from_gcm is set to FALSE !! !! DESCRIPTION : Computes the surface drag coefficient, for cases !! in which it is NOT provided by the coupled atmospheric model LMDZ. The module first uses the !! meteorolgical input to calculate the Richardson Number, which is an indicator of atmospheric !! stability in the surface layer. The formulation used to find this surface drag coefficient is !! dependent on the stability determined. !! !! Designation of wind speed !! !! Calculation of geopotential. This is the definition of Geopotential height (e.g. Jacobson !! eqn.4.47, 2005). (required for calculation of the Richardson Number) !! !! !! Calculation of the virtual air temperature at the surface (required for calculation !! of the Richardson Number) !! !! Calculation of the virtual surface temperature (required for calculation of th !! Richardson Number) !! !! Calculation of the squared wind shear (required for calculation of the Richardson !! Number) !! !! Calculation of the Richardson Number. The Richardson Number is defined as the ratio !! of potential to kinetic energy, or, in the context of atmospheric science, of the !! generation of energy by wind shear against consumption !! by static stability and is an indicator of flow stability (i.e. for when laminar flow !! becomes turbulent and vise versa). It is approximated using the expression below: !! !! The Richardson Number hence calculated is subject to a minimum value: !! !! Computing the drag coefficient. We add the add the height of the vegetation to the !! level height to take into account that the level 'seen' by the vegetation is actually !! the top of the vegetation. Then we we can subtract the displacement height. !! !! For the stable case (i.e $R_i$ $\geq$ 0) !! !! For the unstable case (i.e. $R_i$ < 0) !! !! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere. !! To prevent this, a minimum limit to the drag coefficient is defined as: !! !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): q_cdrag !! !! REFERENCE(S) : !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General !! Circulation Model. Journal of Climate, 6, pp.248-273 !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation !! sur le cycle de l'eau, PhD Thesis, available from: !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf !! - Jacobson M.Z., Fundamentals of Atmospheric Modeling (2nd Edition), published Cambridge !! University Press, ISBN 0-521-54865-9 !! !!_ ================================================================================================================================ SUBROUTINE diffuco_aero (KIDIA, KFDIA, KLON, LDLAND, YDAGS, YDAGF, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, & qsurf, qair, q_cdrag) !! 0. Variable and parameter declaration !! 0.1 Input variables INTEGER(KIND=JPIM), INTENT(IN) :: KIDIA INTEGER(KIND=JPIM), INTENT(IN) :: KFDIA INTEGER(KIND=JPIM), INTENT(IN) :: KLON !! Domain size (-) LOGICAL , INTENT(IN) :: LDLAND(:) TYPE(TAGS) , INTENT(IN) :: YDAGS TYPE(TAGF) , INTENT(IN) :: YDAGF REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1}) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1}) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: zlev !! Height of first atmospheric layer (m) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: z0h !! Surface roughness Length for heat (m) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: z0m !! Surface roughness Length for momentum (m) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: roughheight !! Effective roughness height (m) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: temp_sol !! Ground temperature (K) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: temp_air !! Lowest level temperature (K) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: qsurf !! near surface specific air humidity (kg kg^{-1}) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1}) !! 0.2 Output variables REAL(KIND=JPRB),DIMENSION (KLON), INTENT (out) :: q_cdrag !! Surface drag coefficient (-) !! 0.3 Modified variables !! 0.4 Local variables INTEGER(KIND=JPIM) :: JL, jv REAL(KIND=JPRB) :: speed, zg, zdphi, ztvd, ztvs, zdu2 REAL(KIND=JPRB) :: zri, cd_neut, zscf, cd_tmp REAL(KIND=JPRB) :: wind(KLON) ! Wind (m s^{-1}) REAL(KIND=JPRB) :: cp_h2o ! Specific heat of water vapor (J.kg^{-1}.K^{-1}) REAL(KIND=JPRB) :: rvtmp2 ! Ratio between specific heat of water vapor and dry air ! minus 1 (unitless) REAL(KIND=JPRB) :: retv ! Ratio between molecular weight of dry air and water ! vapor minus 1 (unitless) !_ ================================================================================================================================ !! 1. Initialisation ASSOCIATE(zero=>YDAGF%zero, one=>YDAGF%one, two=>YDAGF%two, three=>YDAGF%three, four=>YDAGF%four, msmlr_h2o=>YDAGF%msmlr_h2o, & & min_wind=>YDAGF%min_wind, cte_grav=>YDAGF%cte_grav, cb=>YDAGF%cb, cc=>YDAGF%cc, cd=>YDAGF%cd, cp_air=>YDAGF%cp_air, & & cepdu2=>YDAGF%cepdu2, ct_karman=>YDAGF%ct_karman, min_qc=>YDAGF%min_qc, & & RMAIR=>YDAGS%RMAIR) cp_h2o = cp_air*(four*YDAGS%RMAIR)/( 3.5_JPRB*msmlr_h2o) rvtmp2 = cp_h2o/cp_air-one retv = RMAIR/msmlr_h2o-one wind(:) = SQRT (u(:)*u(:) + v(:)*v(:)) ! test if we have to work with q_cdrag or to calcul it DO JL=KIDIA,KFDIA !! 1a).1 Designation of wind speed speed = wind(JL) !! 1a).2 Calculation of geopotentiel !! This is the definition of Geopotential height (e.g. Jacobson eqn.4.47, 2005). (required !! for calculation of the Richardson Number) zg = zlev(JL) * cte_grav zdphi = zg/cp_air !! 1a).3 Calculation of the virtual air temperature at the surface !! required for calculation of the Richardson Number ztvd = (temp_air(JL) + zdphi / (one + rvtmp2 * qair(JL))) * (one + retv * qair(JL)) !! 1a).4 Calculation of the virtual surface temperature !! required for calculation of the Richardson Number ztvs = temp_sol(JL) * (one + retv * qsurf(JL)) !! 1a).5 Calculation of the squared wind shear !! required for calculation of the Richardson Number zdu2 = MAX(cepdu2,speed**two) !! 1a).6 Calculation of the Richardson Number !! The Richardson Number is defined as the ratio of potential to kinetic energy, or, in the !! context of atmospheric science, of the generation of energy by wind shear against consumption !! by static stability and is an indicator of flow stability (i.e. for when laminar flow !! becomes turbulent and vise versa). !! It is approximated using the expression below: zri = zg * (ztvd - ztvs) / (zdu2 * ztvd) !! The Richardson Number hence calculated is subject to a minimum value: zri = MAX(MIN(zri,5.),-5.) !! 1a).7 Computing the drag coefficient !! We add the add the height of the vegetation to the level height to take into account !! that the level 'seen' by the vegetation is actually the top of the vegetation. Then we !! we can subtract the displacement height. !! 7.0 Snow smoothering !! Snow induces low levels of turbulence. !! Sensible heat fluxes can therefore be reduced of ~1/3. Pomeroy et al., 1998 cd_neut = ct_karman ** two / ( LOG( (zlev(JL) + roughheight(JL)) / z0m(JL) ) * LOG( (zlev(JL) + roughheight(JL)) / z0h(JL) ) ) !! 1a).7.1 - for the stable case (i.e $R_i$ $\geq$ 0) IF (zri .GE. zero) THEN zscf = SQRT(one + cd * ABS(zri)) cd_tmp=cd_neut/(one + three * cb * zri * zscf) ELSE !! 1a).7.2 - for the unstable case (i.e. $R_i$ < 0) zscf = one / (one + three * cb * cc * cd_neut * SQRT(ABS(zri) * & & ((zlev(JL) + roughheight(JL)) / z0m(JL)))) cd_tmp=cd_neut * (one - three * cb * zri * zscf) ENDIF !! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere. !! To prevent this, a minimum limit to the drag coefficient is defined as: !! q_cdrag(JL) = MAX(cd_tmp, min_qc/MAX(speed,min_wind)) ! In some situations it might be useful to give an upper limit on the cdrag as well. ! The line here should then be uncommented. !q_cdrag(ji) = MIN(q_cdrag(ji), 0.5/MAX(speed,min_wind)) END DO END ASSOCIATE END SUBROUTINE diffuco_aero !! ================================================================================================================================ !! SUBROUTINE : diffuco_inter !! !>\BRIEF This routine computes the partial beta coefficient !! for the interception for each type of vegetation !! !! DESCRIPTION : We first calculate the dry and wet parts of each PFT (wet part = qsintveg/qsintmax). !! The former is submitted to transpiration only (vbeta3 coefficient, calculated in !! diffuco_trans_co2), while the latter is first submitted to interception loss !! (vbeta2 coefficient) and then to transpiration once all the intercepted water has been evaporated !! (vbeta23 coefficient). Interception loss is also submitted to a predictor-corrector test, !! as for snow sublimation. !! !! Calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water !! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss. !! !! !! Calculation of $\beta_3$, the canopy transpiration resistance !! !! We here determine the limitation of interception loss by the water stored on the leaf. !! A first approximation of interception loss is obtained using the old values of !! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc') !! !! !! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction !! 'zqsvegrap'. !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): ::vbeta2, ::vbeta23 !! !! REFERENCE(S) : !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General !! Circulation Model. Journal of Climate, 6, pp. 248-273 !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation !! sur le cycle de l'eau, PhD Thesis, available from: !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf !! - Perrier, A, 1975. Etude physique de l'évaporation dans les conditions naturelles. Annales !! Agronomiques, 26(1-18): pp. 105-123, pp. 229-243 !_ ================================================================================================================================ SUBROUTINE diffuco_inter (KIDIA, KFDIA, KLON, LDLAND, YDAGF, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, & & qsintveg, qsintmax, rstruct, vbeta2, vbeta23) !! 0 Variable and parameter declaration !! 0.1 Input variables INTEGER(KIND=JPIM), INTENT(IN) :: KIDIA INTEGER(KIND=JPIM), INTENT(IN) :: KFDIA INTEGER(KIND=JPIM), INTENT(IN) :: KLON !! Domain size (-) LOGICAL , INTENT(IN) :: LDLAND(:) TYPE(TAGF) , INTENT(IN) :: YDAGF REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1}) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1}) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: rau !! Air Density (kg m^{-3}) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1}) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1}) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: veget !! vegetation fraction for each type (fraction) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2}) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2}) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (in) :: rstruct !! architectural resistance (s m^{-1}) !! 0.2 Output variables REAL(KIND=JPRB),DIMENSION (KLON), INTENT (out) :: vbeta2 !! Beta for interception loss (-) REAL(KIND=JPRB),DIMENSION (KLON), INTENT (out) :: vbeta23 !! Beta for fraction of wetted foliage that will !! transpire (-) !! 0.4 Local variables INTEGER(KIND=JPIM) :: JL, jv !! (-), (-) REAL(KIND=JPRB) :: zqsvegrap, ziltest, zrapp, speed !! REAL(KIND=JPRB) :: wind(KLON) ! Wind (m s^{-1}) !_ ================================================================================================================================ !! 1. Initialize ASSOCIATE(zero=>YDAGF%zero, one=>YDAGF%one, min_wind=>YDAGF%min_wind) vbeta2(:) = zero vbeta23(:) = zero wind(:) = SQRT (u(:)*u(:) + v(:)*v(:)) !! 2. The beta coefficient for interception by vegetation. !DO jv = 2,KVTYPES DO JL=KIDIA,KFDIA IF (veget(JL) .GT. EPSILON(1._JPRM) .AND. qsintveg(JL) .GT. zero ) THEN zqsvegrap = zero IF (qsintmax(JL) .GT. EPSILON(1._JPRM) ) THEN !! !! We calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water !! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss. !! zqsvegrap = MAX(zero, qsintveg(JL) / qsintmax(JL)) END IF speed = MAX(min_wind, wind(JL)) !! Calculation of $\beta_3$, the canopy transpiration resistance vbeta2(JL) = veget(JL) * zqsvegrap * (one / (one + speed * q_cdrag(JL) * rstruct(JL))) !! We here determine the limitation of interception loss by the water stored on the leaf. !! A first approximation of interception loss is obtained using the old values of !! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc') !ziltest = dt_sechiba * vbeta2(JL,jv) * speed * q_cdrag(JL) * rau(JL) * & ziltest = vbeta2(JL) * speed * q_cdrag(JL) * rau(JL) * & & ( qsatt(JL) - qair(JL) ) IF ( ziltest .GT. EPSILON(1._JPRM) ) THEN zrapp = qsintveg(JL) / ziltest IF ( zrapp .LT. one ) THEN !! !! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction !! 'zqsvegrap'. IF ( humrel(JL) >= EPSILON(1._JPRM) ) THEN vbeta23(JL) = MAX(vbeta2(JL) - vbeta2(JL) * zrapp, zero) ELSE ! We don't want transpiration when the soil cannot deliver it vbeta23(JL) = zero ENDIF vbeta2(JL) = vbeta2(JL) * zrapp ENDIF ENDIF END IF ! ! Autre formulation possible pour l'evaporation permettant une transpiration sur tout le feuillage ! !commenter si formulation Nathalie sinon Tristan ! speed = MAX(min_wind, wind(JL)) ! ! vbeta23(JL,jv) = MAX(zero, veget(JL,jv) * (one / (one + speed * q_cdrag(JL) * rstruct(JL,jv))) - vbeta2(JL,jv)) END DO !END DO END ASSOCIATE END SUBROUTINE diffuco_inter ! -------------------------------------------------------------------------- !! SUBROUTINE : diffuco_raerod !! !>\BRIEF Computes the aerodynamic resistance, for cases in which the !! surface drag coefficient is provided by the coupled atmospheric model LMDZ and when the flag !! 'ldq_cdrag_from_gcm' is set to TRUE !! !! DESCRIPTION : Simply computes the aerodynamic resistance, for cases in which the !! surface drag coefficient is provided by the coupled atmospheric model LMDZ. If the surface drag coefficient !! is not provided by the LMDZ or signalled by the flag 'ldq_cdrag_from_gcm' set to FALSE, then the subroutine !! diffuco_aero is called instead of this one. !! !! Calculation of the aerodynamic resistance, for diganostic purposes. First calculate wind speed: ! PARAMETER DESCRIPTION UNITS ! --------- ----------- ----- ! INPUT PARAMETERS (REAL): ! OUTPUT PARAMETERS (REAL): ! *raero* ! REFERENCES ! ---------- ! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations ! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General ! Circulation Model. Journal of Climate, 6, pp.248-273 ! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influence de l'irrigation ! sur le cycle de l'eau, PhD Thesis, available from: ! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf ! -------------------------------------------------------------------------- SUBROUTINE diffuco_raerod (KIDIA, KFDIA, KLON, YDAGF, u, v, q_cdrag, raero) IMPLICIT NONE INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA INTEGER(KIND=JPIM),INTENT(IN) :: KLON TYPE(TAGF) ,INTENT(IN) :: YDAGF REAL(KIND=JPRB) ,INTENT(IN) :: u(:) ! Eastward Lowest level wind velocity (m s^{-1}) REAL(KIND=JPRB) ,INTENT(IN) :: v(:) ! Northward Lowest level wind velocity (m s^{-1}) REAL(KIND=JPRB) ,INTENT(IN) :: q_cdrag(:) ! Surface drag coefficient (-) REAL(KIND=JPRB) ,INTENT(OUT) :: raero(:) ! Aerodynamic resistance (s m^{-1}) !* 0. LOCAL VARIABLES. ! ----- ---------- INTEGER(KIND=JPIM) :: JL ! (-) REAL(KIND=JPRB) :: speed ! (m s^{-1}) REAL(KIND=JPRB) :: wind(KLON) ! Wind (m s^{-1}) ! -------------------------------------------------------------------------- !* 1. Simple calculation of the aerodynamic resistance, for diagnostic purposes ! ------------------------------------------------------------------------- ASSOCIATE(one=>YDAGF%one, min_wind=>YDAGF%min_wind, qsfrict=>YDAGF%qsfrict) wind(:) = SQRT (u(:)*u(:) + v(:)*v(:)) DO JL=KIDIA,KFDIA speed = MAX(min_wind, wind(JL)) raero(JL) = one / (q_cdrag(JL)*speed) ENDDO END ASSOCIATE END SUBROUTINE diffuco_raerod ! -------------------------------------------------------------------------- ! FUNCTION : Arrhenius ! -------------------------------------------------------------------------- FUNCTION Arrhenius (KLON,temp,ref_temp,energy_act,YDAGF) RESULT ( val_arrhenius ) INTEGER(KIND=JPIM),INTENT(IN) :: KLON ! Domain size (-) REAL(KIND=JPRB), INTENT(IN) :: temp(:) ! Temperature (K) REAL(KIND=JPRB), INTENT(IN) :: ref_temp ! Temperature of reference (K) REAL(KIND=JPRB), INTENT(IN) :: energy_act ! Activation Energy (J mol-1) TYPE(TAGF) , INTENT(IN) :: YDAGF !* 0. RESULT ! ------ REAL(KIND=JPRB), DIMENSION(KLON) :: val_arrhenius ! Temperature dependance based on an Arrhenius function (-) val_arrhenius(:)=EXP(((temp(:)-ref_temp)*energy_act)/(ref_temp*YDAGF%RR*(temp(:)))) END FUNCTION Arrhenius ! -------------------------------------------------------------------------- ! FUNCTION : Arrhenius_modified_1d ! -------------------------------------------------------------------------- FUNCTION Arrhenius_modified_1d (KLON,temp,ref_temp,energy_act,energy_deact,entropy,YDAGF) RESULT (val_arrhenius) INTEGER(KIND=JPIM),INTENT(IN) :: KLON ! Domain size (-) REAL(KIND=JPRB), INTENT(IN) :: temp(:) ! Temperature (K) REAL(KIND=JPRB), INTENT(IN) :: ref_temp ! Temperature of reference (K) REAL(KIND=JPRB), INTENT(IN) :: energy_act ! Activation Energy (J mol-1) REAL(KIND=JPRB), INTENT(IN) :: energy_deact ! Deactivation Energy (J mol-1) REAL(KIND=JPRB), INTENT(IN) :: entropy(:) ! Entropy term (J K-1 mol-1) TYPE(TAGF) , INTENT(IN) :: YDAGF !* 0. RESULT ! ------ REAL(KIND=JPRB), DIMENSION(KLON) :: val_arrhenius ! Temperature dependance based on an Arrhenius function (-) val_arrhenius(:)=EXP(((temp(:)-ref_temp)*energy_act)/(ref_temp*YDAGF%RR*(temp(:)))) & * (YDAGF%one + EXP( (ref_temp * entropy(:) - energy_deact) / (ref_temp * YDAGF%RR ))) & / (YDAGF%one + EXP( (temp(:) * entropy(:) - energy_deact) / ( YDAGF%RR*temp(:)))) END FUNCTION Arrhenius_modified_1d ! -------------------------------------------------------------------------- ! FUNCTION : Arrhenius_modified_0d ! -------------------------------------------------------------------------- FUNCTION Arrhenius_modified_0d (KLON, temp,ref_temp,energy_act,energy_deact,entropy,YDAGF) RESULT (val_arrhenius) INTEGER(KIND=JPIM),INTENT(IN) :: KLON ! Domain size (-) REAL(KIND=JPRB), INTENT(IN) :: temp(:) ! Temperature (K) REAL(KIND=JPRB), INTENT(IN) :: ref_temp ! Temperature of reference (K) REAL(KIND=JPRB), INTENT(IN) :: energy_act ! Activation Energy (J mol-1) REAL(KIND=JPRB), INTENT(IN) :: energy_deact ! Deactivation Energy (J mol-1) REAL(KIND=JPRB), INTENT(IN) :: entropy ! Entropy term (J K-1 mol-1) TYPE(TAGF) , INTENT(IN) :: YDAGF !* 0. RESULT ! ------ REAL(KIND=JPRB), DIMENSION(KLON) :: val_arrhenius ! Temperature dependance based on an Arrhenius function (-) val_arrhenius(:)=EXP(((temp(:)-ref_temp)*energy_act)/(ref_temp*YDAGF%RR*(temp(:)))) & * (YDAGF%one + EXP( (ref_temp * entropy - energy_deact) / (ref_temp * YDAGF%RR ))) & / (YDAGF%one + EXP( (temp(:) * entropy - energy_deact) / ( YDAGF%RR*temp(:)))) END FUNCTION Arrhenius_modified_0d !! ================================================================================================================================ !! SUBROUTINE : qsatcalc !! !! This routine calculates the saturated humidity using the pressure !! and the temperature for all pixels. !! !! DESCRIPTION : This routine interpolates qsat between temperatures by the following formula : !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S) : qsat_out !! !! REFERENCE(S) : None !_ ================================================================================================================================ SUBROUTINE qsatcalc (KIDIA,KFDIA,KLON,YDAGF,temp_in,pres_in,qsat_out) IMPLICIT NONE !! 0. Variables and parameters declaration !! 0.1 Input variables INTEGER(KIND=JPIM), INTENT(IN) :: KIDIA INTEGER(KIND=JPIM), INTENT(IN) :: KFDIA INTEGER(KIND=JPIM), INTENT(IN) :: KLON ! Domain size (-) REAL(KIND=JPRB), INTENT(IN) :: temp_in(:) ! Temperature in degre Kelvin (K) REAL(KIND=JPRB), INTENT(IN) :: pres_in(:) ! Pressure (hPa) TYPE(TAGF) , INTENT(IN) :: YDAGF !! 0.2 Output variables REAL(KIND=JPRB), INTENT(OUT) :: qsat_out(:) ! Saturated humidity at the surface (kg of water/kg of air) !! 0.4 Local variables INTEGER(KIND=JPIM) :: jt(KLON) ! Temporary array stocking the truncated temperatures in Kelvin ! (converted into integers) INTEGER(KIND=JPIM) :: ji ! indices (unitless) REAL(KIND=JPRB) :: zz_a(KLON), zz_b(KLON), zz_f(KLON)! Temporary variables INTEGER(KIND=JPIM) :: nbad ! Number of points where the temperature is too high or too low INTEGER(KIND=JPIM),DIMENSION(1):: lo ! Temporary vector to mark the position of the highest temperature ! or the lowest temperature over all the pixels in jt (unitless) INTEGER(KIND=JPIM),PARAMETER :: min_temp=100_JPIM ! Minimum temperature for saturated humidity (K) !_ ================================================================================================================================ !- !! Computes qsat interpolation between two successive temperatures !- jt = INT(temp_in(:)) !! 1. Diagnostic pixels where the temperature is too high nbad = COUNT(jt(:) >= YDAGF%max_temp-1) IF (nbad > 0) THEN WRITE(*,*) ' qsatcalc: temperature too high at ', & & nbad, ' points.' !- !IF (.NOT.diag_qsat) THEN WRITE(*,*) ('WARNING : qsatcalc in FARQUARAHR_MOD, temperature incorect') !ELSE lo = MAXLOC(temp_in(:)) WRITE(*,*) & & 'Maximum temperature ( ',MAXVAL(temp_in),') found at ',lo(1) WHERE (jt(:) >= YDAGF%max_temp-1) jt(:) = YDAGF%max_temp-1 !ENDIF !(.NOT.diag_qsat) !- ENDIF ! (nbad > 0) !! 2. Diagnostic pixels where the temperature is too low nbad = COUNT(jt(:) <= min_temp) IF (nbad > 0) THEN WRITE(*,*) ' qsatcalc: temperature too low at ', & & nbad, ' points.' !- !IF (.NOT.diag_qsat) THEN WRITE(*,*) ('WARNING : qsatcalc in FARQUARAHR_MOD, temperature incorect') !ELSE lo = MINLOC(temp_in(:)) WRITE(*,*) & & 'Minimum temperature ( ',MINVAL(temp_in),') found at ',lo(1) WHERE (jt(:) <= min_temp) jt(:) = min_temp !ENDIF !(.NOT.diag_qsat) !- ENDIF! (nbad > 0) !! 3. Temporary variables needed for interpolation DO ji = KIDIA, KFDIA ! Loop over # pixels zz_f(ji) = temp_in(ji)-FLOAT(jt(ji)) zz_a(ji) = YDAGF%qsfrict(jt(ji)) zz_b(ji) = YDAGF%qsfrict(jt(ji)+1) ENDDO ! Loop over # pixels !- !! 4. Interpolation between these two values !- DO ji = KIDIA, KFDIA ! Loop over # pixels qsat_out(ji) = ((zz_b(ji)-zz_a(ji))*zz_f(ji)+zz_a(ji))/pres_in(ji) ENDDO ! Loop over # pixels END SUBROUTINE qsatcalc END MODULE FARQUHAR_MOD