| 1 | ! ================================================================================================================================= |
|---|
| 2 | ! MODULE : stomate_wind |
|---|
| 3 | ! |
|---|
| 4 | ! CONTACT : yiyingchen@gate.sinica.edu.tw |
|---|
| 5 | ! |
|---|
| 6 | ! LICENCE : IPSL (2006). This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
|---|
| 7 | ! |
|---|
| 8 | !>\BRIEF : Calculates critical wind speed for tree damage for each half-hourly period of the simulation. |
|---|
| 9 | ! This wind speed can be compared to the actual wind speed to decide whether a tree damage occurs. |
|---|
| 10 | !! |
|---|
| 11 | !!\streamlining_n DESCRIPTION: None |
|---|
| 12 | !! |
|---|
| 13 | !! RECENT CHANGE(S): None |
|---|
| 14 | !! |
|---|
| 15 | !! SVN : |
|---|
| 16 | !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/stomate_windthrow.f90 $ |
|---|
| 17 | !! $Date: 2015 -01-22 13:41:13 +0100 (mer., 22 janv. 2014) $ |
|---|
| 18 | !! $Revision: 1831 $ |
|---|
| 19 | !! \streamlining_n |
|---|
| 20 | !_ ================================================================================================================================ |
|---|
| 21 | |
|---|
| 22 | MODULE stomate_windthrow |
|---|
| 23 | |
|---|
| 24 | ! modules used: |
|---|
| 25 | USE ioipsl_para |
|---|
| 26 | USE grid |
|---|
| 27 | USE xios_orchidee |
|---|
| 28 | USE stomate_data |
|---|
| 29 | USE constantes |
|---|
| 30 | USE constantes_soil |
|---|
| 31 | USE pft_parameters |
|---|
| 32 | USE pft_parameters_var |
|---|
| 33 | USE function_library, ONLY: wood_to_dia, Nmax, wood_to_qmdia,wood_to_height, & |
|---|
| 34 | wood_to_cv, cc_to_lai, cc_to_biomass, biomass_to_lai |
|---|
| 35 | |
|---|
| 36 | |
|---|
| 37 | IMPLICIT NONE |
|---|
| 38 | |
|---|
| 39 | ! private & public routines |
|---|
| 40 | |
|---|
| 41 | PRIVATE |
|---|
| 42 | PUBLIC wind_damage |
|---|
| 43 | !!$ wind_clear, streamlining, calculate_force, calculate_bending_moment, & |
|---|
| 44 | !!$ calculate_wind_process, elevate |
|---|
| 45 | |
|---|
| 46 | LOGICAL, SAVE :: firstcall_windthrow = .TRUE. !! First call |
|---|
| 47 | !$OMP THREADPRIVATE(firstcall_windthrow) |
|---|
| 48 | INTEGER(i_std), SAVE :: printlev_loc !! Local level of text output for current module |
|---|
| 49 | !$OMP THREADPRIVATE(printlev_loc) |
|---|
| 50 | |
|---|
| 51 | CONTAINS |
|---|
| 52 | |
|---|
| 53 | |
|---|
| 54 | !! ================================================================================================================================ |
|---|
| 55 | !! SUBROUTINE : wind_clear |
|---|
| 56 | !! |
|---|
| 57 | !>\BRIEF Set the firstcall flag to .TRUE. and activate initialization |
|---|
| 58 | !! |
|---|
| 59 | !_ ================================================================================================================================ |
|---|
| 60 | |
|---|
| 61 | SUBROUTINE wind_clear |
|---|
| 62 | firstcall_windthrow = .TRUE. |
|---|
| 63 | END SUBROUTINE wind_clear |
|---|
| 64 | |
|---|
| 65 | |
|---|
| 66 | !! ================================================================================================================================ |
|---|
| 67 | !! SUBROUTINE : wind_damage |
|---|
| 68 | !! |
|---|
| 69 | !>\BRIEF This subroutine calculates the critical wind speeds for both uprooting and stem breakage |
|---|
| 70 | !! for each half-hourly period of the simulation. The algorithm follows Barry Gardiner's wind |
|---|
| 71 | !! damage risk model called GALES (Hale et al. 2015). |
|---|
| 72 | !! |
|---|
| 73 | !! DESCRIPTION : The purpose of this module is to describe the actual sensitivity of a given PFT to the wind that is present |
|---|
| 74 | !! at the same pixel. For each circumference class in each PFT, two critical wind speed values (CWS) are calculated in every |
|---|
| 75 | !! half-hourly time step by this subroutine, following the approach of differentiating two types of wind damage of trees: |
|---|
| 76 | !! - uprooting (the whole tree is removed from the soil together with its root plate) |
|---|
| 77 | !! - stem breakage (the root plate and the bottom part of the stem remains attached to the ground). |
|---|
| 78 | !! Whichever threshold (i.e. CWS) is reached first for a tree, the according damage type occurs in the model. |
|---|
| 79 | !! |
|---|
| 80 | !! The module completes the following tasks: |
|---|
| 81 | !! 1. Calculates... |
|---|
| 82 | !! |
|---|
| 83 | !! RECENT CHANGE(S): None |
|---|
| 84 | !! |
|---|
| 85 | !! MAIN OUTPUT VARIABLE(S): :: Critical wind speed for stem breakage 10 m above the zero-plane displacement (m/s): U_Break_10 |
|---|
| 86 | !! Critical wind speed for tree overturn 10 m above the zero-plane displacement (m/s): U_Overturn_10 |
|---|
| 87 | !! |
|---|
| 88 | !! REFERENCES : |
|---|
| 89 | !! - Hale, S., Gardiner, B., Nicoll, B., Taylor, P., and Pizzirani, S. (2015). Comparison and validation of three versions |
|---|
| 90 | !! of a forest wind risk model. Environmental Modelling and Software. In Press. |
|---|
| 91 | !! ... |
|---|
| 92 | !! |
|---|
| 93 | !! FLOWCHART : None |
|---|
| 94 | !! \streamlining_n |
|---|
| 95 | !_ ================================================================================================================================ |
|---|
| 96 | |
|---|
| 97 | SUBROUTINE wind_damage (npts, nlevels_tot, circ_class_biomass, & |
|---|
| 98 | veget_max, circ_class_n, wind_speed_daily, plant_status, & |
|---|
| 99 | circ_class_kill, gap_area_save, forest_managed, & |
|---|
| 100 | soil_temp_daily, root_profile, max_wind_speed_storm, count_storm,& |
|---|
| 101 | is_storm) |
|---|
| 102 | |
|---|
| 103 | |
|---|
| 104 | IMPLICIT NONE |
|---|
| 105 | |
|---|
| 106 | !! 0. Variable declarations |
|---|
| 107 | |
|---|
| 108 | !! 0.1 Input variables |
|---|
| 109 | |
|---|
| 110 | INTEGER(i_std), INTENT(in) :: npts !! Domain size @tex $(unitless)$ @endtex |
|---|
| 111 | INTEGER(i_std), INTENT(in) :: nlevels_tot !! Total level numbers in photosynthesis |
|---|
| 112 | INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: forest_managed !! forest management flag (is the forest |
|---|
| 113 | !! being managed?) |
|---|
| 114 | REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass in of an individual tree in each circ class |
|---|
| 115 | !! @tex $(m^{-2})$ @endtex |
|---|
| 116 | REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of trees within each circumference |
|---|
| 117 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT on the ground |
|---|
| 118 | !! May sum to less than unity if the pixel has nobio |
|---|
| 119 | !! area. (unitless, 0-1) |
|---|
| 120 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: plant_status !! Growth and phenological status of the plant |
|---|
| 121 | |
|---|
| 122 | REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: gap_area_save !! Gap area created by more than 30% of the basal area loss |
|---|
| 123 | !! in the last 5 years |
|---|
| 124 | !! Dimension(npts,nvm,wind_years) |
|---|
| 125 | REAL(r_std), DIMENSION(:), INTENT(in) :: soil_temp_daily !! Daily maximum soil temperature at 0.8 meter below ground (K) |
|---|
| 126 | REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: root_profile !! Normalized root mass/length fraction in each soil layer |
|---|
| 127 | !! (0-1, unitless) |
|---|
| 128 | |
|---|
| 129 | !! 0.2 Output variables |
|---|
| 130 | |
|---|
| 131 | !! 0.3 Modified variables |
|---|
| 132 | |
|---|
| 133 | REAL(r_std), DIMENSION(:,:,:,:,:), & |
|---|
| 134 | INTENT(inout) :: circ_class_kill !! Number of trees within a circ that needs |
|---|
| 135 | !! to be killed @tex $(ind m^{-2})$ @endtex |
|---|
| 136 | REAL(r_std), DIMENSION(:), INTENT(inout) :: wind_speed_daily !! Daily maximum wind speed at 2 meter (ms-1) |
|---|
| 137 | |
|---|
| 138 | REAL(r_std), DIMENSION(:), INTENT(inout) :: max_wind_speed_storm !! Daily maximum wind speed at 2 meter (ms-1) during a storm event |
|---|
| 139 | |
|---|
| 140 | INTEGER(i_std), DIMENSION(:), INTENT(inout) :: count_storm !! number of day after a storm |
|---|
| 141 | |
|---|
| 142 | LOGICAL, DIMENSION(:), INTENT(inout) :: is_storm !! are we in a storm event ? |
|---|
| 143 | |
|---|
| 144 | !! 0.4 Local variables |
|---|
| 145 | CHARACTER(30) :: var_name !! To store variable names for I/O |
|---|
| 146 | INTEGER(i_std) :: islm,ipts,ivm !! Index for the loops |
|---|
| 147 | INTEGER(i_std) :: icir,i,iele, iyear !! Index for the loops |
|---|
| 148 | INTEGER(i_std), DIMENSION(npts) :: soil_type !! Soil types: |
|---|
| 149 | !! 1. Free-draining mineral soils; |
|---|
| 150 | !! 2. Gleyed mineral soils; |
|---|
| 151 | !! 3. Peaty mineral soils; and |
|---|
| 152 | !! 4. Deep Peats |
|---|
| 153 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: virtual_circ_class_n !! Virtual number of trees within each circumference |
|---|
| 154 | REAL(r_std), DIMENSION(ncirc,60) :: section_height !! The height of the bottom of the 1-m sections |
|---|
| 155 | !! of the mean tree @tex $(m)$ @endtex. 60 is the |
|---|
| 156 | !! max number of sections thus max tree height is 60m |
|---|
| 157 | REAL(r_std), DIMENSION(ncirc,60) :: branch_width !! Width of canopy) of the 1-m sections of the mean |
|---|
| 158 | !! tree @tex $(m)$ @endtex 60 is the max number of |
|---|
| 159 | !! sections thus max tree height is 60m |
|---|
| 160 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: mean_dbh !! Diameter at breast height (dbh) of the mean tree |
|---|
| 161 | !! of the circumference-class @tex $(m)$ @endtex |
|---|
| 162 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: mean_height !! Height of the mean tree of the given dbh-class |
|---|
| 163 | !! @tex $(m)$ @endtex |
|---|
| 164 | REAL(r_std), DIMENSION(ncirc) :: lai !! Leaf area index @tex $(unitless)$ @endtex |
|---|
| 165 | REAL(r_std), DIMENSION(ncirc) :: canopy_breadth !! Canopy breadth (width of canopy) @tex $(m)$ @endtex |
|---|
| 166 | REAL(r_std), DIMENSION(ncirc) :: canopy_depth !! Canopy depth (height of canopy) @tex $(m)$ @endtex |
|---|
| 167 | REAL(r_std), DIMENSION(ncirc) :: canopy_height !! The height of the starting point of the canopy from |
|---|
| 168 | !! below @tex $(m)$ @endtex |
|---|
| 169 | REAL(r_std), DIMENSION(ncirc) :: mid_canopy !! The middle height of the canopy @tex $(m)$ @endtex |
|---|
| 170 | REAL(r_std), DIMENSION(ncirc) :: mean_dbh_height !! The height of the diameter (by default it is at |
|---|
| 171 | !! breast height = 1.3 m) @tex $(m)$ @endtex |
|---|
| 172 | REAL(r_std), DIMENSION(ncirc) :: snow_mass !! Weight of the snow on the tree @tex $(kg)$ @endtex |
|---|
| 173 | REAL(r_std), DIMENSION(ncirc) :: diameter_c2 !! Variable for calculating diameters |
|---|
| 174 | !! @tex $(unitless)$ @endtex |
|---|
| 175 | REAL(r_std), DIMENSION(ncirc) :: diameter_canopy_base !! Diameter at the base of the canopy |
|---|
| 176 | !! @tex $(m)$ @endtex |
|---|
| 177 | REAL(r_std), DIMENSION(ncirc) :: diameter_c1 !! Variable for calculating diameters |
|---|
| 178 | !! @tex $(unitless)$ @endtex |
|---|
| 179 | REAL(r_std), DIMENSION(ncirc) :: crown_depth !! Crown depth of the tree in each circumference class |
|---|
| 180 | REAL(r_std), DIMENSION(ncirc) :: crown_width !! Crown width of the mean tree (m) |
|---|
| 181 | REAL(r_std), DIMENSION(ncirc) :: crown_volume !! Crown volume of the mean tree @tex $(m^{3})$ @endtex |
|---|
| 182 | REAL(r_std), DIMENSION(ncirc) :: crown_mass !! Mass of the mean tree crown @tex $(kg)$ @endtex |
|---|
| 183 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: stem_mass !! Mass of the mean tree stem @tex $(kg)$ @endtex |
|---|
| 184 | REAL(r_std), DIMENSION(ncirc) :: current_spacing !! Spacing between trees in the given dbh-class |
|---|
| 185 | !! @tex $(m)$ @endtexi |
|---|
| 186 | INTEGER(i_std), DIMENSION(npts,nvm) :: rooting_depth !! Rooting depth of the mean tree (three categories; |
|---|
| 187 | !! 1: shallow, 2: deep, 3:average |
|---|
| 188 | !! See constantes_soil_var.f90 for the definitions |
|---|
| 189 | REAL(r_std) :: streamlining_c !! Streamlining parameter. @tex $(unitless)$ @endtex. |
|---|
| 190 | !! Streamlining is the change of shape of the crowns |
|---|
| 191 | !! due to wind. |
|---|
| 192 | REAL(r_std) :: streamlining_n !! Streamlining parameter. @tex $(unitless)$ @endtex. |
|---|
| 193 | !! Streamlining is the change of shape of the crowns |
|---|
| 194 | !! due to wind. |
|---|
| 195 | REAL(r_std) :: area_timber_removals_5_years !! Area of the total timber removal in the |
|---|
| 196 | !! previous 5 year @tex $(m^{2})$ @endtex |
|---|
| 197 | REAL(r_std), DIMENSION(ncirc) :: tree_heights_from_edge !! Width of the forest edge area, i.e. area_closer (in number of tree heights) |
|---|
| 198 | !! @tex $(unitless)$ @endtex |
|---|
| 199 | REAL(r_std) :: n_gap !! The number of gaps with an area of clear_cut_max @tex $(unitless)$ @endtex |
|---|
| 200 | REAL(r_std) :: length_side_gap !! The length of one side of the square-shaped gap @tex $(m)$ @endtex |
|---|
| 201 | REAL(r_std), DIMENSION(ncirc) :: area_around_gap !! The framing edge area around the gap with different gustiness |
|---|
| 202 | !! @tex $(m^{2})$ @endtex |
|---|
| 203 | REAL(r_std), DIMENSION(ncirc) :: area_total_closer !! The total area of forest edge areas affected by wind @tex $(m^{2})$ @endtex |
|---|
| 204 | REAL(r_std), DIMENSION(ncirc) :: area_total_further !! The total area of forest without edge areas @tex $(m^{2})$ @endtex |
|---|
| 205 | REAL(r_std), DIMENSION(ncirc) :: ratio_spacing_height !! Ratio of tree height and tree spacing @tex $(unitless)$ @endtex |
|---|
| 206 | REAL(r_std) :: gap_size !! The length (width of the cross-section) of the gap @tex $(m)$ @endtex |
|---|
| 207 | REAL(r_std), DIMENSION(ncirc) :: mean_gap_factor !! Variable used for calculating gustiness |
|---|
| 208 | !! @tex $(unitless)$ @endtex |
|---|
| 209 | REAL(r_std), DIMENSION(ncirc) :: max_gap_factor !! Variable used for calculating gustiness |
|---|
| 210 | !! @tex $(unitless)$ @endtex |
|---|
| 211 | REAL(r_std), DIMENSION(ncirc) :: term_a !! Variable used for calculating gustiness |
|---|
| 212 | !! @tex $(unitless)$ @endtex |
|---|
| 213 | REAL(r_std), DIMENSION(ncirc) :: term_b !! Variable used for calculating gustiness |
|---|
| 214 | !! @tex $(unitless)$ @endtex |
|---|
| 215 | REAL(r_std), DIMENSION(ncirc) :: old_gust_closer !! Variable used for calculating gustiness |
|---|
| 216 | !! @tex $(unitless)$ @endtex |
|---|
| 217 | REAL(r_std), DIMENSION(ncirc) :: old_gust_further !! Variable used for calculating gustiness |
|---|
| 218 | !! @tex $(unitless)$ @endtex |
|---|
| 219 | REAL(r_std), DIMENSION(ncirc) :: old_gust_edge !! Variable used for calculating gustiness |
|---|
| 220 | !! @tex $(unitless)$ @endtex |
|---|
| 221 | REAL(r_std), DIMENSION(ncirc) :: old_gust_edge_gap !! Variable used for calculating gustiness |
|---|
| 222 | !! @tex $(unitless)$ @endtex |
|---|
| 223 | REAL(r_std), DIMENSION(ncirc) :: new_gust_closer !! Variable used for calculating gustiness |
|---|
| 224 | !! @tex $(unitless)$ @endtex |
|---|
| 225 | REAL(r_std), DIMENSION(ncirc) :: new_gust_further !! Variable used for calculating gustiness |
|---|
| 226 | !! @tex $(unitless)$ @endtex |
|---|
| 227 | REAL(r_std), DIMENSION(ncirc) :: new_gust_edge !! Variable used for calculating gustiness |
|---|
| 228 | !! @tex $(unitless)$ @endtex |
|---|
| 229 | REAL(r_std), DIMENSION(ncirc) :: new_gust_edge_gap_closer !! Variable used for calculating gustiness |
|---|
| 230 | !! @tex $(unitless)$ @endtex |
|---|
| 231 | REAL(r_std), DIMENSION(ncirc) :: new_gust_edge_gap_further !! Variable used for calculating gustiness |
|---|
| 232 | !! @tex $(unitless)$ @endtex |
|---|
| 233 | REAL(r_std), DIMENSION(ncirc) :: mean_bm_gust_factor_closer !! Mean bending moment of the forest edge area |
|---|
| 234 | !! (for calculating the edge factor) |
|---|
| 235 | !! @tex $(Nm/kg)$ @endtex |
|---|
| 236 | REAL(r_std), DIMENSION(ncirc) :: mean_bm_gust_factor_further !! Mean bending moment of the inner forest area (for calculating |
|---|
| 237 | !! the edge factor) @tex $(Nm/kg)$ @endtex |
|---|
| 238 | REAL(r_std), DIMENSION(ncirc) :: edge_factor !! The ratio of the mean bending moment in the forest edge area |
|---|
| 239 | !! and mean bending moment in the inner forest area |
|---|
| 240 | !! @tex $(unitless)$ @endtex |
|---|
| 241 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: gust_factor_closer !! Gust factor for calculating the critical wind speed of the forest |
|---|
| 242 | !! edge area @tex $(unitless)$ @endtex |
|---|
| 243 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: gust_factor_further !! Gust factor for calculating the critical wind speed of the inner |
|---|
| 244 | !! forest area @tex $(unitless)$ @endtex |
|---|
| 245 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: g_closer !! Gustiness of the forest edge area @tex $(unitless)$ @endtex |
|---|
| 246 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: g_further !! Gustiness of the inner forest area @tex $(unitless)$ @endtex |
|---|
| 247 | REAL(r_std), DIMENSION(npts,nvm) :: overturning_moment_multiplier!! Overturning moment multiplier representing the effects of species, |
|---|
| 248 | !! soil type, soil depths and senescence @tex $(Nm/kg)$ @endtex |
|---|
| 249 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: max_overturning_moment !! The maximum overturning moment the mean tree can withstand |
|---|
| 250 | !! @tex $(Nm/kg)$ @endtex |
|---|
| 251 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: max_breaking_moment !! The maximum breaking moment the mean tree can withstand |
|---|
| 252 | !! @tex $(Nm/kg)$ @endtex |
|---|
| 253 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: overturning_moment_closer !! The maximum overturning moment the mean tree can withstand in the |
|---|
| 254 | !! forest edge area with gustiness taken into account @tex $(Nm/kg)$ @endtex |
|---|
| 255 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: overturning_moment_further !! The maximum overturning moment the mean tree can withstand in the |
|---|
| 256 | !! inner forest area with gustiness taken into account @tex $(Nm/kg)$ @endtex |
|---|
| 257 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: breaking_moment_closer !! The maximum breaking moment the mean tree can withstand in the |
|---|
| 258 | !! forest edge area with gustiness taken into account @tex $(Nm/kg)$ @endtex |
|---|
| 259 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: breaking_moment_further !! The maximum breaking moment the mean tree can withstand in the |
|---|
| 260 | !! inner forest area with gustiness taken into account @tex $(Nm/kg)$ @endtex |
|---|
| 261 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_break_closer !! Critical wind speed for stem breakage in the forest edge area |
|---|
| 262 | !! @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 263 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_break_closer_10 !! Critical wind speed for stem breakage in the forest edge area |
|---|
| 264 | !! @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 265 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_overturn_closer !! Critical wind speed for overturning in the forest edge area |
|---|
| 266 | !! @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 267 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_overturn_closer_10 !! Critical wind speed for overturning in the forest edge area |
|---|
| 268 | !! @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 269 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_final_closer !! The critical wind speed of the final form of damage (breakage / uprooting) |
|---|
| 270 | !! in the forest edge area @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 271 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: d_h_ratio !! Array for storing the stand spacing and tree height ratio for output (-) |
|---|
| 272 | REAL(r_std), DIMENSION(npts) :: wind_speed_actual !! The actual wind speed over the pixel |
|---|
| 273 | !! @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 274 | LOGICAL :: wind_damage_type_uprooting_closer !! The final form of damage (breakage / uprooting) in the forest edge area |
|---|
| 275 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: wind_damage_rate_closer !! Damage rate (%) in the forest edge area in case of damage |
|---|
| 276 | !! @tex $(0-1)$ @endtex |
|---|
| 277 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_break_further !! Critical wind speed for stem breakage in the inner forest area |
|---|
| 278 | !! @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 279 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_break_further_10 !! Critical wind speed for stem breakage in the inner forest area |
|---|
| 280 | !! @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 281 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_overturn_further !! Critical wind speed for overturning in the inner forest area |
|---|
| 282 | !! @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 283 | !Jina: Seems it if for 10 m above the z0 but has the same description. |
|---|
| 284 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_overturn_further_10 !! Critical wind speed for overturning in the inner forest area |
|---|
| 285 | !! @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 286 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_final_further !! The critical wind speed of the final form of damage (breakage / uprooting) |
|---|
| 287 | !! in the inner forest area @tex $(m/s; mean hourly wind speed)$ @endtex |
|---|
| 288 | LOGICAL :: wind_damage_type_uprooting_further !! The final form of damage (breakage / uprooting) in the inner forest area |
|---|
| 289 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: wind_damage_rate_further !! Damage rate (%) in the forest edge area in case |
|---|
| 290 | !! of damage in the inner forest area @tex $(0-1)$ @endtex |
|---|
| 291 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: z0_wind !! surface roughness calculated from the subroutine streamlining() |
|---|
| 292 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: zhd_wind !! zero plane displacement height calculated from the subroutine streamlining(m) |
|---|
| 293 | INTEGER(i_std), DIMENSION(npts,nvm,ncirc) :: wind_damage_type_closer !! Integer store the wind damage type for closer area ibreakage=1, ioverturning=2 |
|---|
| 294 | INTEGER(i_std), DIMENSION(npts,nvm,ncirc) :: wind_damage_type_further !! Integer store the wind damage type for further area |
|---|
| 295 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: kill_break !! Biomass loss through stem breakage from storm |
|---|
| 296 | !! @tex $(ind m^{-2})$ @endtex |
|---|
| 297 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: kill_uproot !! Number of trees killed through uproot from storm |
|---|
| 298 | !! @tex $(ind m^{-2})$ @endtex |
|---|
| 299 | REAL(r_std) :: crown_top !! Temporal variable for storing the tallest tree crown height in different PFTs |
|---|
| 300 | REAL(r_std) :: crown_bottom !! Temporal variable for storing the smallest tree crown height in different PFTs |
|---|
| 301 | REAL(r_std) :: distance_from_force !! The distance from the force was set to 1.3 for breakage and 0.0 for overturning |
|---|
| 302 | REAL(r_std), DIMENSION(ncirc) :: closer_area_fraction !! Closer area fraction |
|---|
| 303 | REAL(r_std), DIMENSION(ncirc) :: further_area_fraction !! Further area fraction |
|---|
| 304 | REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: biomass |
|---|
| 305 | |
|---|
| 306 | ! Temp output Variables for GALES model intercomparision |
|---|
| 307 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: canopy_depth_out |
|---|
| 308 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: tree_h_edge_out |
|---|
| 309 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: mean_gap_f_out |
|---|
| 310 | REAL(r_std), DIMENSION(npts,nvm) :: gap_size_out |
|---|
| 311 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: edge_factor_out |
|---|
| 312 | REAL(r_std) :: roots_below !! temporary variable containing the share of roots below 80 cm depth (0-1, unitless) |
|---|
| 313 | !_ ============================================================================================================================================================================ |
|---|
| 314 | |
|---|
| 315 | IF (firstcall_windthrow) THEN |
|---|
| 316 | !! Initialize local printlev |
|---|
| 317 | printlev_loc=get_printlev('stomate_windthrow') |
|---|
| 318 | |
|---|
| 319 | ! The code is specific for CRU-NCEP. If a different forcing |
|---|
| 320 | ! is used, the Max Wind Ratio parameters should be adjusted |
|---|
| 321 | CALL ipslerr_p(2,'Wind throw module assumes that CRU-NCEP forcing is used',& |
|---|
| 322 | '','','') |
|---|
| 323 | END IF |
|---|
| 324 | |
|---|
| 325 | IF (printlev_loc>=2) WRITE(numout,*) 'Entering stomate_windthrow.f90' |
|---|
| 326 | |
|---|
| 327 | !! 0. Initialized the values of variables |
|---|
| 328 | wind_damage_rate_closer(:,:,:) = zero |
|---|
| 329 | wind_damage_rate_further(:,:,:) = zero |
|---|
| 330 | cws_final_further(:,:,:) = zero |
|---|
| 331 | cws_final_closer(:,:,:) = zero |
|---|
| 332 | d_h_ratio(:,:,:) = zero |
|---|
| 333 | kill_uproot = zero |
|---|
| 334 | kill_break = zero |
|---|
| 335 | cws_break_further_10 = zero |
|---|
| 336 | cws_break_closer_10 = zero |
|---|
| 337 | cws_overturn_further_10 = zero |
|---|
| 338 | cws_overturn_closer_10 = zero |
|---|
| 339 | virtual_circ_class_n = zero |
|---|
| 340 | |
|---|
| 341 | !-Debug- |
|---|
| 342 | IF (printlev_loc .GT. 4) THEN |
|---|
| 343 | DO ipts = 1, npts |
|---|
| 344 | WRITE(numout,*) 'Daily maximum windspeed passed from stomate to stomate_winthrow:', wind_speed_daily(ipts) |
|---|
| 345 | WRITE(numout,*) 'Daily maximum soil temperature at 80cm below ground:', soil_temp_daily(ipts) |
|---|
| 346 | ENDDO |
|---|
| 347 | ENDIF |
|---|
| 348 | !- |
|---|
| 349 | |
|---|
| 350 | !! 1. Getting the input variables in the required form for the actual subroutine |
|---|
| 351 | ! calculate the total biomass |
|---|
| 352 | biomass(:,:,:,:) = cc_to_biomass(npts,nvm,circ_class_biomass,circ_class_n) |
|---|
| 353 | |
|---|
| 354 | !! 1.1 Convert 6-hr CRU-NCEP into half-hourly wind speed |
|---|
| 355 | ! The max_wind_ratio parameters were determined by relating fluxnet |
|---|
| 356 | ! half-hourly data to 6h CRU-NCEP data. For this reason the |
|---|
| 357 | ! parameters should be considered specific to the CRU-NCEP forcing. |
|---|
| 358 | ! The max_wind_ratio parameters were also found to depend on the |
|---|
| 359 | ! spatial aggregation, hence, they are specific to the 0.5 degree |
|---|
| 360 | ! CRU-NCEP data (see Chen et al for details). |
|---|
| 361 | |
|---|
| 362 | ! The following relationship to convert CRU-NCEP into half-hourly |
|---|
| 363 | ! wind speeds was found: |
|---|
| 364 | !wind_speed_daily(:) = -5.922*wind_speed_daily(:) & |
|---|
| 365 | ! + 2.051*wind_speed_daily(:)**2. & |
|---|
| 366 | ! - 0.191*wind_speed_daily(:)**3. & |
|---|
| 367 | ! + 0.006*wind_speed_daily(:)**4. |
|---|
| 368 | !WHERE (wind_speed_daily(:) .LE. zero) |
|---|
| 369 | ! wind_speed_daily(:) = zero |
|---|
| 370 | !ENDWHERE |
|---|
| 371 | |
|---|
| 372 | ! In case this convertion did not result in the observed wind |
|---|
| 373 | ! damage, we kept a linear tuning parameter. Ideally, no tuning |
|---|
| 374 | ! should be applied and therefore daily_max_tune should be 1.0 |
|---|
| 375 | !wind_speed_daily(:) = wind_speed_daily(:) * daily_max_tune |
|---|
| 376 | |
|---|
| 377 | ! Storms do not last for exactly 24 hours. |
|---|
| 378 | ! Storms that pass overnight would be accounted for twice |
|---|
| 379 | ! unless we add criteria to tell ORCHIDEE this is just a |
|---|
| 380 | ! single storm. The value wind_speed_damage_thr |
|---|
| 381 | ! (set to 20.0 m/s) corresponds to the wind speed |
|---|
| 382 | ! threshold above which is_storm flag is set the TRUE |
|---|
| 383 | ! and nb_days_storm (set to 5.0 days) correspond to the |
|---|
| 384 | ! number of days at which the max wind speed is less than the |
|---|
| 385 | ! threshold which means that the storm is over and hence |
|---|
| 386 | ! is_storm flag to FALSE. As a consequence the windstorm |
|---|
| 387 | ! damage is calculated nb_days_storm (set to 5 days) |
|---|
| 388 | ! after the end of a storm. One main consequences of this |
|---|
| 389 | ! new approach is that damages due to wind is only accounted |
|---|
| 390 | ! when the flag is_storm is set to TRUE. |
|---|
| 391 | WHERE (wind_speed_daily(:) .GT. & |
|---|
| 392 | wind_speed_storm_thr .AND. .NOT. is_storm(:)) |
|---|
| 393 | is_storm(:) = .TRUE. |
|---|
| 394 | count_storm(:) = 0 |
|---|
| 395 | ENDWHERE |
|---|
| 396 | |
|---|
| 397 | DO ipts = 1, npts |
|---|
| 398 | |
|---|
| 399 | if(is_storm(ipts)) then |
|---|
| 400 | max_wind_speed_storm(ipts)= & |
|---|
| 401 | max(max_wind_speed_storm(ipts), & |
|---|
| 402 | wind_speed_daily(ipts)) |
|---|
| 403 | if(wind_speed_daily(ipts) .LT. & |
|---|
| 404 | wind_speed_storm_thr) then |
|---|
| 405 | count_storm(ipts) = count_storm(ipts) + 1 |
|---|
| 406 | endif |
|---|
| 407 | if(count_storm(ipts) .LT. nb_days_storm) then |
|---|
| 408 | cycle |
|---|
| 409 | endif |
|---|
| 410 | count_storm(ipts) = 0 |
|---|
| 411 | is_storm(:) = .FALSE. |
|---|
| 412 | else |
|---|
| 413 | cycle |
|---|
| 414 | endif |
|---|
| 415 | |
|---|
| 416 | !! 1.3 Defining the soil type |
|---|
| 417 | ! This sub-routine uses 4 types of soil for calculating the |
|---|
| 418 | ! critical wind speed, so 12 different soil types present in |
|---|
| 419 | ! other modules of ORCHIDEE need to be grouped into these |
|---|
| 420 | ! 4 generic types. |
|---|
| 421 | |
|---|
| 422 | ! +++CHECK+++ |
|---|
| 423 | ! For the moment only 3 soil types are distinguished. |
|---|
| 424 | ! Needs to be checked with all 12 USDA soil types. |
|---|
| 425 | ! If more soil type are used it will become possible |
|---|
| 426 | ! to make use of the different soils distinguished in |
|---|
| 427 | ! the windthrow module. The code could then look as |
|---|
| 428 | ! follow: |
|---|
| 429 | !!$ SELECT CASE() |
|---|
| 430 | !!$ CASE (x:y) |
|---|
| 431 | !!$ soil_type(ipts) = ifreedrainage |
|---|
| 432 | !!$ CASE (x:y) |
|---|
| 433 | !!$ soil_type(ipts) = igleyed |
|---|
| 434 | !!$ CASE (x:y) |
|---|
| 435 | !!$ soil_type(ipts) = ipeaty |
|---|
| 436 | !!$ CASE (x:y) |
|---|
| 437 | !!$ soil_type(ipts) = ipeat |
|---|
| 438 | !!$ CASE DEFAULT |
|---|
| 439 | !!$ CALL ipslerr_p(3, 'stomate_windthrow.f90',& |
|---|
| 440 | !!$ 'soil type is not defined','','') |
|---|
| 441 | !!$ END SELECT |
|---|
| 442 | |
|---|
| 443 | ! For the moment all soils are assumed to be freely draining |
|---|
| 444 | soil_type(ipts) = ifree_draining |
|---|
| 445 | ! +++++++++++++ |
|---|
| 446 | |
|---|
| 447 | DO ivm = 2, nvm |
|---|
| 448 | |
|---|
| 449 | ! Initilaize |
|---|
| 450 | area_timber_removals_5_years = zero |
|---|
| 451 | n_gap = zero |
|---|
| 452 | |
|---|
| 453 | ! Don't do the calculations if there is no vegetation |
|---|
| 454 | ! or the vegetation is bare soil, grass or crop |
|---|
| 455 | IF ((veget_max(ipts,ivm).LT.min_stomate) .OR. .NOT.is_tree(ivm)) THEN |
|---|
| 456 | CYCLE |
|---|
| 457 | ENDIF |
|---|
| 458 | |
|---|
| 459 | !! 1.4 Calculating the rooting depth |
|---|
| 460 | ! If root density (proportion of root mass) is higher |
|---|
| 461 | ! than 40% below 80 cm below the soil surface, then we |
|---|
| 462 | ! treat rooting depth as "deep". Otherwise, it's "shallow". |
|---|
| 463 | ! Note that when using the static root profile deep vs shallow |
|---|
| 464 | ! will be the same for all pixels belonhing to the same PFT |
|---|
| 465 | ! beacuse it solely depends on humcste. The threshold of 40% |
|---|
| 466 | ! was chosen to ensure that all trees are shallow rooting with |
|---|
| 467 | ! the static root profile. |
|---|
| 468 | IF (grnd_80.EQ.nslm) THEN |
|---|
| 469 | ! Handle the case where the deepest soil layer is less |
|---|
| 470 | ! than 80 cm. All roots must then be above 80 cm. |
|---|
| 471 | roots_below = un |
|---|
| 472 | ELSE |
|---|
| 473 | ! Calculate the share of roots below 80 cm |
|---|
| 474 | roots_below = zero |
|---|
| 475 | DO islm = grnd_80,nslm |
|---|
| 476 | roots_below = roots_below + root_profile(ipts,ivm,islm,istruc) |
|---|
| 477 | END DO |
|---|
| 478 | END IF |
|---|
| 479 | IF (roots_below > 0.40) THEN |
|---|
| 480 | rooting_depth(ipts,ivm) = ideep |
|---|
| 481 | ELSE |
|---|
| 482 | rooting_depth(ipts,ivm) = ishallow |
|---|
| 483 | END IF |
|---|
| 484 | |
|---|
| 485 | !--Debug |
|---|
| 486 | IF (printlev_loc >= 4) THEN |
|---|
| 487 | WRITE(numout,*) 'rooting depth ', rooting_depth(ipts,ivm) |
|---|
| 488 | ENDIF |
|---|
| 489 | |
|---|
| 490 | !! 1.5 Set crown parameters |
|---|
| 491 | ! To set crown parameters leaf mass needs to be available |
|---|
| 492 | IF (SUM(circ_class_biomass(ipts,ivm,:,ileaf,icarbon)) .GT. min_stomate) THEN |
|---|
| 493 | ! There is a canopy |
|---|
| 494 | streamlining_c = streamlining_c_leaf(ivm) |
|---|
| 495 | streamlining_n = streamlining_n_leaf(ivm) |
|---|
| 496 | ELSE |
|---|
| 497 | ! If the tree is leafless, the following parameters |
|---|
| 498 | ! become different. The weight of leaves will be absent. |
|---|
| 499 | ! Streamlining changes as the surface of leaves is absent. |
|---|
| 500 | streamlining_c = streamlining_c_leafless(ivm) |
|---|
| 501 | streamlining_n = streamlining_n_leafless(ivm) |
|---|
| 502 | END IF |
|---|
| 503 | |
|---|
| 504 | !-Debug- |
|---|
| 505 | IF(printlev_loc .GT. 4 .AND. ivm .EQ. test_pft)THEN |
|---|
| 506 | WRITE(numout,*) 'Initialisation of critical windspeed' |
|---|
| 507 | WRITE(numout,*) 'ipts, ivm, ',ipts, ivm |
|---|
| 508 | WRITE(numout,*) 'soil_type(ipts), ',soil_type(ipts) |
|---|
| 509 | WRITE(numout,*) 'rooting_depth_type,(1:deep, 2:shallow,3:average), ',& |
|---|
| 510 | rooting_depth(ipts,ivm) |
|---|
| 511 | WRITE(numout,*) 'streamlining_c, ', streamlining_c |
|---|
| 512 | WRITE(numout,*) 'streamlining_n, ', streamlining_n |
|---|
| 513 | ENDIF |
|---|
| 514 | !- |
|---|
| 515 | |
|---|
| 516 | !! 1.6 Calculate stand characteristics |
|---|
| 517 | ! Calculating diameter at breast height (dbh) in m of the |
|---|
| 518 | ! mean tree of all circ classes. |
|---|
| 519 | mean_dbh(ipts,ivm,:) = wood_to_dia(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm) |
|---|
| 520 | |
|---|
| 521 | ! Calculating the height in m of the mean tree of the given |
|---|
| 522 | ! dbh-class using |
|---|
| 523 | mean_height(ipts,ivm,:) = wood_to_height(& |
|---|
| 524 | circ_class_biomass(ipts,ivm,:,:,icarbon),ivm) |
|---|
| 525 | |
|---|
| 526 | ! Crown_volume is calculated as the woody biomass of an individual tree |
|---|
| 527 | crown_volume(:) = wood_to_cv( circ_class_biomass(ipts,ivm,:,:,icarbon), & |
|---|
| 528 | circ_class_n(ipts,ivm,:),ivm) |
|---|
| 529 | |
|---|
| 530 | ! Crown_width is calculated as(crown_vloume * 6.0/pi)^(1/3) |
|---|
| 531 | crown_width(:) = ((crown_volume(:)*6.0)/pi)**(1.0/3.0) |
|---|
| 532 | |
|---|
| 533 | ! ORCHIDEE assumes that the crowns are spherical. Hence, width |
|---|
| 534 | ! breath and depth are all the same. Set canopy_breadth as the |
|---|
| 535 | ! crown_width for all circumference classes |
|---|
| 536 | canopy_breadth(:) = crown_width(:) |
|---|
| 537 | canopy_depth(:) = crown_width(:) |
|---|
| 538 | canopy_height(:) = mean_height(ipts,ivm,:) |
|---|
| 539 | |
|---|
| 540 | ! Calculate gaps and the surrounding areas |
|---|
| 541 | ! IMPORTANT: AS WE DO NOT HAVE INFORMATION ON STAND EDGES IN THE |
|---|
| 542 | ! ORCHIDEE SIMULATIONS, A POSSIBLE PROXY TO USE HERE IS A VARIABLE |
|---|
| 543 | ! THAT REPRESENTS EARLIER |
|---|
| 544 | |
|---|
| 545 | ! GAP AREAS IN THE PREVIOUS 5 YEARS. AFTER |
|---|
| 546 | ! SUCH A PERIOD, THE FOREST EDGE IS ASSUMED TO GET ACCLIMATIZED |
|---|
| 547 | ! TO THE NEW CIRCUMSTANCES, AND IT'S VULNERABILITY DIMINISHES |
|---|
| 548 | ! (see Persson, 1975; Valinger and Fridman, 2011) |
|---|
| 549 | ! Note that gap area for the current year is not calculated yet (which is |
|---|
| 550 | ! calculated at the end of the year in stomate_season), so the gap |
|---|
| 551 | ! created t-1 to t-5 is applied. |
|---|
| 552 | ! Calculating the total area of timber removals in the previous |
|---|
| 553 | ! five years (set by the variable ::wind_years |
|---|
| 554 | area_timber_removals_5_years = zero |
|---|
| 555 | DO iyear = 1, wind_years |
|---|
| 556 | area_timber_removals_5_years = area_timber_removals_5_years + & |
|---|
| 557 | gap_area_save(ipts,ivm,iyear) |
|---|
| 558 | END DO |
|---|
| 559 | |
|---|
| 560 | ! Number of gaps with the area of the maximum allowed clearfelling: |
|---|
| 561 | n_gap = area_timber_removals_5_years / clear_cut_max |
|---|
| 562 | |
|---|
| 563 | ! The length of the sides of the squared-shaped gaps: |
|---|
| 564 | length_side_gap = SQRT(clear_cut_max) |
|---|
| 565 | |
|---|
| 566 | !-Debug- |
|---|
| 567 | IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN |
|---|
| 568 | WRITE(numout,*) 'crown_volumn, ', crown_volume(:) |
|---|
| 569 | WRITE(numout,*) 'crown_width, ', crown_width(:) |
|---|
| 570 | WRITE(numout,*) 'crown_depth, ', canopy_depth(:) |
|---|
| 571 | WRITE(numout,*) 'canopy_height, ', mean_height(:,:,:) |
|---|
| 572 | ENDIF |
|---|
| 573 | !- |
|---|
| 574 | |
|---|
| 575 | ! loop over the circ class |
|---|
| 576 | DO icir = 1, ncirc |
|---|
| 577 | |
|---|
| 578 | !-Debug- |
|---|
| 579 | IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN |
|---|
| 580 | WRITE(numout,*) 'biomass, isapabove, ', biomass(ipts,ivm,isapabove,icarbon) |
|---|
| 581 | WRITE(numout,*) 'biomass, iheartabove, ', biomass(ipts,ivm,iheartabove,icarbon) |
|---|
| 582 | WRITE(numout,*) 'circ_class_biomass, isapabove, ', & |
|---|
| 583 | circ_class_biomass(ipts,ivm,icir,isapabove,icarbon) |
|---|
| 584 | WRITE(numout,*) 'circ_class_biomass, iheatabove, ', & |
|---|
| 585 | circ_class_biomass(ipts,ivm,icir,iheartabove,icarbon) |
|---|
| 586 | ENDIF |
|---|
| 587 | !- |
|---|
| 588 | |
|---|
| 589 | !! 1.6 Calculating the current tree spacing in the PFT |
|---|
| 590 | ! The current tree spacing between tree stems are calculated |
|---|
| 591 | ! from the number of stems per hectare. circ_class_n is the |
|---|
| 592 | ! number of stems per m2; it's converted to stems per hectare |
|---|
| 593 | ! with m2_to_ha which has a value of 10000. |
|---|
| 594 | ! The virtual tree spacing assumes that the LAI remains |
|---|
| 595 | ! constant and the calculates how many trees of a given |
|---|
| 596 | ! circumference class would be needed to result in the set LAI. |
|---|
| 597 | ! Because lai changes throughout the season because of leaf |
|---|
| 598 | ! fall, wood mass was used instead. Because LAI is related to |
|---|
| 599 | ! wood mass (see allocation). The principle remains the same |
|---|
| 600 | IF ( (circ_class_biomass(ipts,ivm,icir,isapabove,icarbon) + & |
|---|
| 601 | circ_class_biomass(ipts,ivm,icir,iheartabove,icarbon)) >= zero ) THEN |
|---|
| 602 | virtual_circ_class_n(ipts,ivm,icir) = & |
|---|
| 603 | (biomass(ipts,ivm,isapabove,icarbon) + & |
|---|
| 604 | biomass(ipts,ivm,iheartabove,icarbon)) / & |
|---|
| 605 | (circ_class_biomass(ipts,ivm,icir,isapabove,icarbon) + & |
|---|
| 606 | circ_class_biomass(ipts,ivm,icir,iheartabove,icarbon)) |
|---|
| 607 | ELSE |
|---|
| 608 | virtual_circ_class_n(ipts,ivm,icir) = zero |
|---|
| 609 | ENDIF |
|---|
| 610 | |
|---|
| 611 | ! Now the spacing is calculated based on the stem density |
|---|
| 612 | ! based on the virtual circ class to beeter deal with |
|---|
| 613 | ! unmanaged forest and heterogenous canopy. |
|---|
| 614 | IF (virtual_circ_class_n(ipts,ivm,icir) .GT. min_stomate) THEN |
|---|
| 615 | ! The original equation is current_spacing(icir) = 100.0 / & |
|---|
| 616 | ! SQRT(m2_to_ha * virtual_circ_class_n(ipts,ivm,icir)) |
|---|
| 617 | ! given that m2_to_ha = 10000 it reduces to: |
|---|
| 618 | current_spacing(icir) = 1.0 / SQRT(virtual_circ_class_n(ipts,ivm,icir)) |
|---|
| 619 | ELSE |
|---|
| 620 | ! In principle, we should never enter this condition. The |
|---|
| 621 | ! spacing between 2 trees equals the size of the pixel. It is |
|---|
| 622 | ! thus empty. We defined this case such that if we do end here, |
|---|
| 623 | ! a very large cws will be calculated and the code should not |
|---|
| 624 | ! suffer from divide by zero. |
|---|
| 625 | current_spacing(icir) = SQRT(area(ipts)) |
|---|
| 626 | WRITE(numout,*) 'WARNING from WINDFALL module:' |
|---|
| 627 | WRITE(numout,*) 'the current spacing set to pixel spacing:', current_spacing(icir) |
|---|
| 628 | CALL ipslerr_p(2,'Wind throw module',& |
|---|
| 629 | 'current spacing was set to pixel dimensions',& |
|---|
| 630 | 'it is therefore assumed to be empty','') |
|---|
| 631 | ENDIF |
|---|
| 632 | |
|---|
| 633 | !! 1.7 Calculate green stem mass |
|---|
| 634 | ! Stem mass is calculated as the aboveground woody biomass |
|---|
| 635 | ! minus the branches. circ_class_biomass is the biomass of |
|---|
| 636 | ! an individual tree (gC/tree) |
|---|
| 637 | ! Here, we convert the unit to (kgC/tree) by dividing kilo_to_unit (1000) |
|---|
| 638 | stem_mass(ipts,ivm,icir) = (circ_class_biomass(ipts,ivm,icir,isapabove,icarbon) + & |
|---|
| 639 | circ_class_biomass(ipts,ivm,icir,iheartabove,icarbon)) * & |
|---|
| 640 | (un - branch_ratio(ivm)) / kilo_to_unit |
|---|
| 641 | |
|---|
| 642 | !-Debug- |
|---|
| 643 | IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN |
|---|
| 644 | WRITE(numout,*) 'stem_mass,',stem_mass(ipts,ivm,icir) |
|---|
| 645 | WRITE(numout,*) 'green_density,', green_density(ivm) |
|---|
| 646 | WRITE(numout,*) 'pipe_sensity,', pipe_density(ivm) |
|---|
| 647 | WRITE(numout,*) 'kilo to unit,', kilo_to_unit |
|---|
| 648 | WRITE(numout,*) 'virtual_circ_class', virtual_circ_class_n(ipts,ivm,icir) |
|---|
| 649 | WRITE(numout,*) 'sla,',sla(ivm) |
|---|
| 650 | END IF |
|---|
| 651 | !- |
|---|
| 652 | |
|---|
| 653 | ! Wind throw occurs on living trees so the stem mass should be |
|---|
| 654 | ! converted from dry biomass in gC to wet biomass in kg to be used |
|---|
| 655 | ! in the GALES equations. Stem_mass thus converted as follows: |
|---|
| 656 | ! (green_density / (pipe_density/1000)) from dry wood mass |
|---|
| 657 | ! (kgC/tree) to wet wood mass (Kg/trees) |
|---|
| 658 | stem_mass(ipts,ivm,icir) = stem_mass(ipts, ivm,icir) * & |
|---|
| 659 | green_density(ivm)/(pipe_density(ivm)/kilo_to_unit) |
|---|
| 660 | |
|---|
| 661 | !! 1.8 Calculating leaf area index |
|---|
| 662 | ! We need to use the virtual biomass in each circumference class |
|---|
| 663 | ! to calculate lai, and tree_heights_from_edge was calculated as |
|---|
| 664 | ! function of lai. Use the function biomass_to_lai to account for |
|---|
| 665 | ! a dynamic lai calculation. |
|---|
| 666 | lai(icir) = biomass_to_lai(circ_class_biomass(ipts,ivm,icir,ileaf,icarbon)*& |
|---|
| 667 | virtual_circ_class_n(ipts,ivm,icir),ivm) |
|---|
| 668 | |
|---|
| 669 | !-Debug- |
|---|
| 670 | IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN |
|---|
| 671 | WRITE(numout,*) 'check spacing and circ_number in a diameter class' |
|---|
| 672 | WRITE(numout,*) 'virtual_tree_numbers:', & |
|---|
| 673 | virtual_circ_class_n(ipts,ivm,icir), ' for circ class:', icir |
|---|
| 674 | WRITE(numout,*) 'current spacing: ',& |
|---|
| 675 | current_spacing(icir),'mean diameter: ', mean_dbh(ipts,ivm,icir) |
|---|
| 676 | WRITE(numout,*) 'area:', area(ipts) |
|---|
| 677 | WRITE(numout,*) 'ipts, ',ipts,'ivm, ',ivm,'icir, ',icir |
|---|
| 678 | WRITE(numout,*) 'mean dbh, ', mean_dbh(ipts,ivm,icir) |
|---|
| 679 | WRITE(numout,*) 'mean height, ', mean_height(ipts,ivm,icir) |
|---|
| 680 | WRITE(numout,*) 'stem_mass, ', stem_mass(ipts,ivm,icir) |
|---|
| 681 | WRITE(numout,*) 'crown_volume, ',crown_volume(icir) |
|---|
| 682 | WRITE(numout,*) 'crown_mass, ',crown_mass(icir) |
|---|
| 683 | WRITE(numout,*) 'crown_width, ',crown_width(icir) |
|---|
| 684 | WRITE(numout,*) 'lai, ',lai(icir) |
|---|
| 685 | WRITE(numout,*) 'virtual_circ_n, ', virtual_circ_class_n(ipts,ivm,icir) |
|---|
| 686 | WRITE(numout,*) 'real_circ_n, ', circ_class_n(ipts,ivm,icir) |
|---|
| 687 | WRITE(numout,*) 'virtual_spacing, ',current_spacing(icir) |
|---|
| 688 | WRITE(numout,*) 'real_spacing, ', SQRT(1.0/circ_class_n(ipts,ivm,icir)) |
|---|
| 689 | WRITE(numout,*) 'End of Section 1: Calculated canopy structure & |
|---|
| 690 | & in difffernt circumferences for each PFT' |
|---|
| 691 | ENDIF |
|---|
| 692 | ! |
|---|
| 693 | |
|---|
| 694 | !! 2. Calculating the gustiness of wind in the PFTs |
|---|
| 695 | |
|---|
| 696 | !! 2.1 Calculate gaps and the surrounding areas |
|---|
| 697 | ! IMPORTANT: AS WE DO NOT HAVE INFORMATION ON STAND EDGES IN THE |
|---|
| 698 | ! ORCHIDEE SIMULATIONS, A POSSIBLE PROXY TO USE HERE IS A VARIABLE |
|---|
| 699 | ! THAT REPRESENTS EARLIER HARVESTS IN THE PREVIOUS 5 YEARS. AFTER |
|---|
| 700 | ! SUCH A PERIOD, THE FOREST EDGE IS ASSUMED TO GET ACCLIMATIZED |
|---|
| 701 | ! TO THE NEW CIRCUMSTANCES, AND IT'S VULNERABILITY DIMINISHES |
|---|
| 702 | ! (see Persson, 1975; Valinger and Fridman, 2011) |
|---|
| 703 | ! There is a big increase in gustiness after a certain distance from |
|---|
| 704 | ! the edge of the forest, hence we divide the area of the PFT into two |
|---|
| 705 | ! areas: |
|---|
| 706 | ! - the area which is further than X tree heights from any edge that |
|---|
| 707 | ! was created in the previous five years (area_total_further) |
|---|
| 708 | ! - the area which is closer than X tree heights from any edge that |
|---|
| 709 | ! was created in the previous five years (area_total_closer). |
|---|
| 710 | ! Of course, these smaller areas depend on the fragmentation of the |
|---|
| 711 | ! gaps (created in the last five years). As we have no information |
|---|
| 712 | ! on this, we assume a fixed gap size representing the most likely |
|---|
| 713 | ! size of clear cuts (clear_cut_max) and that the |
|---|
| 714 | ! gaps are square-shaped. Then, we divide the area of harvest in the |
|---|
| 715 | ! previous five years (area_timber_removals_5_years) by clear_cut_max. |
|---|
| 716 | ! Hence, from the number of gaps we can calculate the area_total_closer. |
|---|
| 717 | ! As each circumference class will have its forest area divided, we |
|---|
| 718 | ! calculate two gustiness values and hence, two critical wind speeds. |
|---|
| 719 | ! The value of X depends on the leaf area index of the PFT and is |
|---|
| 720 | ! calculated as follows: |
|---|
| 721 | ! tree_heights_from_edge is defined as "X/h", X sets to 28/LAI and h |
|---|
| 722 | ! is set to meant_tree_height. |
|---|
| 723 | ! This comes from the equation of calculating canopy penetration depth |
|---|
| 724 | ! in Pan et al. 2017?. Ying Pan work on the large eddy simulation and |
|---|
| 725 | ! is now at Pennsilvania State University |
|---|
| 726 | |
|---|
| 727 | ! Here, lai(icir) for each diameter class is calculated by virtual_class_n |
|---|
| 728 | ! a threshold "0.1" is used to avoid the issue of divided by small lai < 0.1 |
|---|
| 729 | ! for the cases of deciduous tree species |
|---|
| 730 | IF (lai(icir) .GT. 0.1) THEN |
|---|
| 731 | tree_heights_from_edge(icir) = (28.0/lai(icir)) / mean_height(ipts,ivm,icir) |
|---|
| 732 | ELSE |
|---|
| 733 | tree_heights_from_edge(icir) = (28.0/0.1) / mean_height(ipts,ivm,icir) |
|---|
| 734 | ENDIF |
|---|
| 735 | |
|---|
| 736 | ! As suggested by Barry (the father of GALES), the tree height from edge |
|---|
| 737 | ! should be limited to 9.0 |
|---|
| 738 | IF (tree_heights_from_edge(icir) .GT. 9.0) tree_heights_from_edge(icir) = 9.0 |
|---|
| 739 | |
|---|
| 740 | ! The frame-shaped area around the clear cut gap with X as the width |
|---|
| 741 | ! of the frame: |
|---|
| 742 | area_around_gap(icir) = ((length_side_gap + & |
|---|
| 743 | tree_heights_from_edge(icir)*mean_height(ipts,ivm,icir)*2.0)**2.0) - & |
|---|
| 744 | clear_cut_max |
|---|
| 745 | |
|---|
| 746 | ! Total area of the "frames" in the particular dbh-class of the particula |
|---|
| 747 | ! PFT:the multiplication by 0.25 is because of the assumption that wind |
|---|
| 748 | ! direction does not change during the half-hourly time step, therefore |
|---|
| 749 | ! only one of the four sides of the gaps are affected. |
|---|
| 750 | ! Use MIN to make sure area_total_closer and area_total_further are |
|---|
| 751 | ! positive. |
|---|
| 752 | ! +++JINA This area should not be below zero. I added zero to MIN |
|---|
| 753 | ! for now, but maybe need to add error to avoid miss |
|---|
| 754 | ! calculations. |
|---|
| 755 | area_total_closer(icir) = MIN(n_gap * area_around_gap(icir) * 0.25, & |
|---|
| 756 | area(ipts)*veget_max(ipts,ivm) - & |
|---|
| 757 | (n_gap * clear_cut_max), zero) |
|---|
| 758 | |
|---|
| 759 | ! The total area which is nor gap neither "frame": |
|---|
| 760 | area_total_further(icir) = area(ipts) * veget_max(ipts,ivm) - & |
|---|
| 761 | area_total_closer(icir) - (n_gap * clear_cut_max) |
|---|
| 762 | |
|---|
| 763 | IF (area(ipts).GT.zero) THEN |
|---|
| 764 | ! Get the area fraction |
|---|
| 765 | closer_area_fraction(icir) = area_total_closer(icir) / (area(ipts)*veget_max(ipts,ivm)) |
|---|
| 766 | further_area_fraction(icir) = area_total_further(icir) / (area(ipts)*veget_max(ipts,ivm)) |
|---|
| 767 | ELSE |
|---|
| 768 | WRITE(numout,*) 'ERROR: the total area of the pixel is zero' |
|---|
| 769 | CALL ipslerr_p(3,'stomate_windthrow',& |
|---|
| 770 | 'area of the pixel is zero','','') |
|---|
| 771 | ENDIF |
|---|
| 772 | |
|---|
| 773 | IF(ok_wlsk .AND. ivm .EQ. test_pft .AND. ipts .EQ. test_grid) THEN |
|---|
| 774 | WRITE(numout,*) 'wlskwlsk_area, ipts, icir',ipts, icir |
|---|
| 775 | WRITE(numout,*) 'n_gap',n_gap |
|---|
| 776 | WRITE(numout,*) 'area',area(ipts) |
|---|
| 777 | WRITE(numout,*) 'area-harvested',area(ipts)-(n_gap*clear_cut_max) |
|---|
| 778 | WRITE(numout,*) 'gap_area',gap_area_save(ipts,ivm,:) |
|---|
| 779 | WRITE(numout,*) 'closer_area_f',closer_area_fraction(icir) |
|---|
| 780 | WRITE(numout,*) 'further_area_f',further_area_fraction(icir) |
|---|
| 781 | ENDIF |
|---|
| 782 | |
|---|
| 783 | !-Debug- |
|---|
| 784 | IF(printlev_loc .GT. 4 .AND. ivm .EQ. test_pft)THEN |
|---|
| 785 | WRITE(numout,*) 'Section 2.0: Calculating the gustiness of wind & |
|---|
| 786 | & in the PFT and area fractions in the grid.' |
|---|
| 787 | WRITE(numout,*) 'area_closer_fraction, ', closer_area_fraction(icir) |
|---|
| 788 | WRITE(numout,*) 'area_further_fraction, ', further_area_fraction(icir) |
|---|
| 789 | ENDIF |
|---|
| 790 | !- |
|---|
| 791 | |
|---|
| 792 | !! 2.2 Calculate ratio to spacing height |
|---|
| 793 | ! The ratio of spacing to height describes the functional (for |
|---|
| 794 | ! wind) stand density. If there are many tall trees, the spacing |
|---|
| 795 | ! is low and thus ratio_spacing_height will be low as well. |
|---|
| 796 | ! Values of the ratio of mean tree spacing and height are limited |
|---|
| 797 | ! between 0.075 and 0.45. Bigger values don't make any difference, |
|---|
| 798 | ! smaller values are not really possible in real life. |
|---|
| 799 | ! The D/h value are presented in Hale et al., 2015 |
|---|
| 800 | ratio_spacing_height(icir) = MAX(MIN(current_spacing(icir) / & |
|---|
| 801 | mean_height(ipts,ivm,icir),0.45),0.075) |
|---|
| 802 | |
|---|
| 803 | ! Copy the ratio_spacing_height to the D_H_ratio for output |
|---|
| 804 | d_h_ratio(ipts,ivm,icir) = ratio_spacing_height(icir) |
|---|
| 805 | |
|---|
| 806 | !! 2.3 Calculate functional gap size |
|---|
| 807 | ! Calculate whether the gap is large enough to increase |
|---|
| 808 | ! the sensitivity of the remaining trees to wind throw. |
|---|
| 809 | ! GALES uses the length of the edge of the square-shaped |
|---|
| 810 | ! gaps calculated in the previous section. |
|---|
| 811 | gap_size = length_side_gap |
|---|
| 812 | |
|---|
| 813 | ! If the gap is wider than 10 times the tree height, |
|---|
| 814 | ! then wind loading on the stand edge does not increase any |
|---|
| 815 | ! more. The threshold of 10 is used in the equations below |
|---|
| 816 | IF (gap_size > 10.0 * mean_height(ipts,ivm,icir)) THEN |
|---|
| 817 | mean_gap_factor(icir) = 1. |
|---|
| 818 | max_gap_factor(icir) = 1. |
|---|
| 819 | ELSE |
|---|
| 820 | mean_gap_factor(icir) = (0.001 * (gap_size / & |
|---|
| 821 | mean_height(ipts,ivm,icir))**0.562) / (0.001 * 10.**0.562) |
|---|
| 822 | max_gap_factor(icir) = (0.0064 * (gap_size / & |
|---|
| 823 | mean_height(ipts,ivm,icir))**0.3467) / (0.0064 * 10.**0.3467) |
|---|
| 824 | END IF |
|---|
| 825 | |
|---|
| 826 | ! In the following section, there are three types of terms for forest edges: |
|---|
| 827 | ! the edge (edge), the area around the edge (closer) and the area further |
|---|
| 828 | ! away from the edge (further). Only closer and further are used in |
|---|
| 829 | ! the calculation of cws. Gustiness is different in both areas. |
|---|
| 830 | |
|---|
| 831 | ! Calculating term_a and term_b for the new gust method. The definition of |
|---|
| 832 | ! the term_a and term_b are in the Hale et al. (2015) eq(7) |
|---|
| 833 | ! D/h is the ratio_spacing_height |
|---|
| 834 | ! x/h is the tree_heights_from_edge |
|---|
| 835 | term_a(icir) = MAX(-2.1 * ratio_spacing_height(icir) + 0.91, zero) |
|---|
| 836 | term_b(icir) = 1.0611 * LOG(ratio_spacing_height(icir)) + 4.2 |
|---|
| 837 | |
|---|
| 838 | ! New gust factors: |
|---|
| 839 | ! The gustiness away from the edge is always higher than at the edge. |
|---|
| 840 | ! We applied term_a*tree_heights_from_edge + term_b for gustiness calculation. |
|---|
| 841 | ! The gustiness away from the edge is always higher than at the edge. |
|---|
| 842 | ! For the further area, tree_height_from_edge is always greater than 9.0 |
|---|
| 843 | ! thus, the above calculation will produce an even higher value for gustiness |
|---|
| 844 | ! calculation. So, We set tree_heights_from_edge as 9.0. |
|---|
| 845 | ! The gustiness at the further area will than be higher then the gustiness at |
|---|
| 846 | ! the closer area. |
|---|
| 847 | new_gust_closer(icir) = term_a(icir) * tree_heights_from_edge(icir) + term_b(icir) |
|---|
| 848 | new_gust_further(icir) = term_a(icir) * 9.0 + term_b(icir) |
|---|
| 849 | new_gust_edge(icir) = term_a(icir) * zero + term_b(icir) |
|---|
| 850 | |
|---|
| 851 | ! Calculating gustiness: |
|---|
| 852 | g_closer(ipts,ivm,icir) = new_gust_closer(icir) |
|---|
| 853 | g_further(ipts,ivm,icir) = new_gust_further(icir) |
|---|
| 854 | |
|---|
| 855 | ! The following section is used to calculate the edge_factor. This section calculates |
|---|
| 856 | ! how the mean loading on a tree changes as a function of tree height, spacing, distance |
|---|
| 857 | ! from edge and size of any upwind gap. An upwind gap greater than 10 tree heights |
|---|
| 858 | ! is assumed to be an infinite gap. The output (edge_factor) is used to modify the |
|---|
| 859 | ! calculated wind loading on the tree, which assumes a steady wind and a position well |
|---|
| 860 | ! inside the forest. |
|---|
| 861 | mean_bm_gust_factor_closer(icir) = (0.68 * ratio_spacing_height(icir) - & |
|---|
| 862 | 0.0385) + (-0.68 * ratio_spacing_height(icir) + 0.4785) * ((1.7239 * & |
|---|
| 863 | ratio_spacing_height(icir) + 0.0316)**tree_heights_from_edge(icir)) * & |
|---|
| 864 | mean_gap_factor(icir) |
|---|
| 865 | mean_bm_gust_factor_further(icir) = (0.68 * ratio_spacing_height(icir) - & |
|---|
| 866 | 0.0385) + (-0.68 * ratio_spacing_height(icir) + 0.4785) * ((1.7239 * & |
|---|
| 867 | ratio_spacing_height(icir) + 0.0316)**9.0) * mean_gap_factor(icir) |
|---|
| 868 | edge_factor(icir) = mean_bm_gust_factor_closer(icir) / mean_bm_gust_factor_further(icir) |
|---|
| 869 | ! Copy array for output |
|---|
| 870 | edge_factor_out(ipts,ivm,icir) = edge_factor(icir) |
|---|
| 871 | ! Calculating gust_factor accounting for gustiness and edge_factor: |
|---|
| 872 | gust_factor_closer(ipts,ivm,icir) = g_closer(ipts,ivm,icir) * edge_factor(icir) |
|---|
| 873 | |
|---|
| 874 | ! Suggested by Barry |
|---|
| 875 | ! For further we force edge_factor = 1. This is because everything is |
|---|
| 876 | ! compared to wind load on tree in the middle of the forest |
|---|
| 877 | edge_factor(icir) = 1.0 |
|---|
| 878 | gust_factor_further(ipts,ivm,icir) = g_further(ipts,ivm,icir) * edge_factor(icir) |
|---|
| 879 | |
|---|
| 880 | !-Debug- |
|---|
| 881 | IF(printlev_loc .GT. 4 .AND. ivm .EQ. test_pft)THEN |
|---|
| 882 | WRITE(numout,*) 'End of Section 2.0: Calculating the & |
|---|
| 883 | & gustiness of wind in the PFT' |
|---|
| 884 | WRITE(numout,*) 'term_a, ', term_a(icir) |
|---|
| 885 | WRITE(numout,*) 'term_b, ', term_b(icir) |
|---|
| 886 | WRITE(numout,*) 'new_gust_closer, ', new_gust_closer(icir) |
|---|
| 887 | WRITE(numout,*) 'ratio_spacing_height, ', ratio_spacing_height(:) |
|---|
| 888 | WRITE(numout,*) 'new_gust_further, ', new_gust_further(icir) |
|---|
| 889 | WRITE(numout,*) 'new_gust_edge, ', new_gust_edge(icir) |
|---|
| 890 | WRITE(numout,*) 'tree_heights_from_edge, ', tree_heights_from_edge(icir) |
|---|
| 891 | WRITE(numout,*) 'Correction factor for the forest & |
|---|
| 892 | & edge dimension in number of circumference' |
|---|
| 893 | WRITE(numout,*) 'edge_factor, ', edge_factor(icir) |
|---|
| 894 | WRITE(numout,*) 'gustness_closer, ', g_closer(ipts,ivm,icir) |
|---|
| 895 | WRITE(numout,*) 'gustness_futher, ', g_further(ipts,ivm,icir) |
|---|
| 896 | WRITE(numout,*) 'gust_factor_closer, ' , gust_factor_closer(ipts,ivm,icir) |
|---|
| 897 | WRITE(numout,*) 'gust_factor_further, ', gust_factor_further(ipts,ivm,icir) |
|---|
| 898 | WRITE(numout,*) 'Edge factor:', edge_factor(icir) |
|---|
| 899 | WRITE(numout,*) 'Tree height from edge:', tree_heights_from_edge(icir) |
|---|
| 900 | WRITE(numout,*) 'Mean bm_gust_factor_closer:', & |
|---|
| 901 | mean_bm_gust_factor_closer(icir) |
|---|
| 902 | WRITE(numout,*) 'Mean bm_gust_factor_further:', & |
|---|
| 903 | mean_bm_gust_factor_further(icir) |
|---|
| 904 | WRITE(numout,*) 'Mean gap_factor:', mean_gap_factor(icir) |
|---|
| 905 | ENDIF |
|---|
| 906 | !- |
|---|
| 907 | |
|---|
| 908 | !! 3. Calculating critical moments (tree mechanics section) |
|---|
| 909 | ! overturning_moment_multiplier: this comes form the species parameters file and |
|---|
| 910 | ! is related to the depth and type of the soil. Rooting depth has three categories: |
|---|
| 911 | ! Shallow (less then 0.8 m deep); Deep (deeper than 0.8 m); Average (average values |
|---|
| 912 | ! to be used when rooting depth is unknown). Soil type has four categories: |
|---|
| 913 | ! Free-draining mineral soils; Gleyed mineral soils; Peaty mineral soils; Deep peats. |
|---|
| 914 | ! As these are generic soil types, a modified soil map of Europe would be needed |
|---|
| 915 | ! for the optimal use of WINDTHROW. |
|---|
| 916 | IF (soil_type(ipts) == ifree_draining .AND. Rooting_depth(ipts,ivm) == ishallow & |
|---|
| 917 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 918 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 919 | overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_shallow(ivm) |
|---|
| 920 | ELSE IF (soil_type(ipts) == ifree_draining .AND. rooting_depth(ipts,ivm) == ishallow & |
|---|
| 921 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 922 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 923 | overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_shallow_leafless(ivm) |
|---|
| 924 | ELSE IF (soil_type(ipts) == ifree_draining .AND. rooting_depth(ipts,ivm) == ideep & |
|---|
| 925 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR.& |
|---|
| 926 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 927 | overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_deep(ivm) |
|---|
| 928 | ELSE IF (soil_type(ipts) == ifree_draining .AND. rooting_depth(ipts,ivm) == ideep & |
|---|
| 929 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 930 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 931 | overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_deep_leafless(ivm) |
|---|
| 932 | ELSE IF (soil_type(ipts) == ifree_draining .AND. rooting_depth(ipts,ivm) == iaverage & |
|---|
| 933 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 934 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 935 | overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_average(ivm) |
|---|
| 936 | ELSE IF (soil_type(ipts) == ifree_draining .AND. rooting_depth(ipts,ivm) == iaverage & |
|---|
| 937 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 938 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 939 | overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_average_leafless(ivm) |
|---|
| 940 | ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == ishallow & |
|---|
| 941 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 942 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 943 | overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_shallow(ivm) |
|---|
| 944 | ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == ishallow & |
|---|
| 945 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 946 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 947 | overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_shallow_leafless(ivm) |
|---|
| 948 | ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == ideep & |
|---|
| 949 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 950 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 951 | overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_deep(ivm) |
|---|
| 952 | ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == ideep & |
|---|
| 953 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 954 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 955 | overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_deep_leafless(ivm) |
|---|
| 956 | ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == iaverage & |
|---|
| 957 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 958 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 959 | overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_average(ivm) |
|---|
| 960 | ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == iaverage & |
|---|
| 961 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 962 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 963 | overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_average_leafless(ivm) |
|---|
| 964 | ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == ishallow & |
|---|
| 965 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 966 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 967 | overturning_moment_multiplier(ipts,ivm) = overturning_peaty_shallow(ivm) |
|---|
| 968 | ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == ishallow & |
|---|
| 969 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 970 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 971 | overturning_moment_multiplier(ipts,ivm) = overturning_peaty_shallow_leafless(ivm) |
|---|
| 972 | ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == ideep & |
|---|
| 973 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 974 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 975 | overturning_moment_multiplier(ipts,ivm) = overturning_peaty_deep(ivm) |
|---|
| 976 | ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == ideep & |
|---|
| 977 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 978 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 979 | overturning_moment_multiplier(ipts,ivm) = overturning_peaty_deep_leafless(ivm) |
|---|
| 980 | ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == iaverage & |
|---|
| 981 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 982 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 983 | overturning_moment_multiplier(ipts,ivm) = overturning_peaty_average(ivm) |
|---|
| 984 | ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == iaverage & |
|---|
| 985 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 986 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 987 | overturning_moment_multiplier(ipts,ivm) = overturning_peaty_average_leafless(ivm) |
|---|
| 988 | ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == ishallow & |
|---|
| 989 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 990 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 991 | overturning_moment_multiplier(ipts,ivm) = overturning_peat_shallow(ivm) |
|---|
| 992 | ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == ishallow & |
|---|
| 993 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 994 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 995 | overturning_moment_multiplier(ipts,ivm) = overturning_peat_shallow_leafless(ivm) |
|---|
| 996 | ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == ideep & |
|---|
| 997 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 998 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 999 | overturning_moment_multiplier(ipts,ivm) = overturning_peat_deep(ivm) |
|---|
| 1000 | ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == ideep & |
|---|
| 1001 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 1002 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 1003 | overturning_moment_multiplier(ipts,ivm) = overturning_peat_deep_leafless(ivm) |
|---|
| 1004 | ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == iaverage & |
|---|
| 1005 | .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & |
|---|
| 1006 | plant_status(ipts,ivm).EQ.ipresenescence)) THEN |
|---|
| 1007 | overturning_moment_multiplier(ipts,ivm) = overturning_peat_average(ivm) |
|---|
| 1008 | ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == iaverage & |
|---|
| 1009 | .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & |
|---|
| 1010 | plant_status(ipts,ivm).NE.ipresenescence)) THEN |
|---|
| 1011 | overturning_moment_multiplier(ipts,ivm) = overturning_peat_average_leafless(ivm) |
|---|
| 1012 | ELSE |
|---|
| 1013 | !average value if parameters are missing |
|---|
| 1014 | overturning_moment_multiplier(ipts,ivm) = 125.0 |
|---|
| 1015 | END IF |
|---|
| 1016 | |
|---|
| 1017 | ! Maximum moment calculated for tree overturning. |
|---|
| 1018 | max_overturning_moment(ipts,ivm,icir) = overturning_moment_multiplier(ipts,ivm) * & |
|---|
| 1019 | stem_mass(ipts,ivm,icir) |
|---|
| 1020 | |
|---|
| 1021 | ! Maximum moment calculated for stem breakage. |
|---|
| 1022 | max_breaking_moment(ipts,ivm,icir) = modulus_rupture(ivm) * f_knot(ivm) * pi * & |
|---|
| 1023 | (mean_dbh(ipts,ivm,icir)**3.0) / 32.0 |
|---|
| 1024 | |
|---|
| 1025 | ! This is where the gustiness of wind is accounted for so that half hourly mean |
|---|
| 1026 | ! wind speeds can be used in the simulations. |
|---|
| 1027 | ! In the end, moments and hence, critical wind speeds are calculated |
|---|
| 1028 | ! for both areas (closer & further) with different gustiness. |
|---|
| 1029 | overturning_moment_closer(ipts,ivm,icir) = max_overturning_moment(ipts,ivm,icir) / & |
|---|
| 1030 | gust_factor_closer(ipts,ivm,icir) |
|---|
| 1031 | overturning_moment_further(ipts,ivm,icir) = max_overturning_moment(ipts,ivm,icir) / & |
|---|
| 1032 | gust_factor_further(ipts,ivm,icir) |
|---|
| 1033 | breaking_moment_closer(ipts,ivm,icir) = max_breaking_moment(ipts,ivm,icir) / & |
|---|
| 1034 | gust_factor_closer(ipts,ivm,icir) |
|---|
| 1035 | breaking_moment_further(ipts,ivm,icir) = max_breaking_moment(ipts,ivm,icir) / & |
|---|
| 1036 | gust_factor_further(ipts,ivm,icir) |
|---|
| 1037 | |
|---|
| 1038 | !-Debug- |
|---|
| 1039 | IF(printlev_loc .GT. 4 .AND. ivm .EQ. test_pft)THEN |
|---|
| 1040 | WRITE(numout,*) 'PFT type, ', ivm |
|---|
| 1041 | WRITE(numout,*) 'POINT id, ', ipts |
|---|
| 1042 | WRITE(numout,*) 'SENSCENCE type,', plant_status(ipts,ivm) |
|---|
| 1043 | WRITE(numout,*) 'Overturning multiplier for free drainage average root leafless all PFTs,', & |
|---|
| 1044 | overturning_free_draining_average_leafless(:) |
|---|
| 1045 | WRITE(numout,*) 'Overturning multiplier for free drainage average root leaf all PFTs,', & |
|---|
| 1046 | overturning_free_draining_average(:) |
|---|
| 1047 | WRITE(numout,*) 'Overturning multiplier for free drainage shallow root leafless all PFTs,', & |
|---|
| 1048 | overturning_free_draining_shallow_leafless(:) |
|---|
| 1049 | WRITE(numout,*) 'Overturning multiplier for free drainage shallow root leaf all PFTs,', & |
|---|
| 1050 | overturning_free_draining_shallow(:) |
|---|
| 1051 | WRITE(numout,*) 'Overturning multiplier for free drainage deep root leafless all PFTs,', & |
|---|
| 1052 | overturning_free_draining_deep_leafless(:) |
|---|
| 1053 | WRITE(numout,*) 'Overturning multiplier for free drainage deep root leaf all PFTs,', & |
|---|
| 1054 | overturning_free_draining_deep(:) |
|---|
| 1055 | WRITE(numout,*) 'End of Section 3 : Calculating critical moments (tree mechanics section)' |
|---|
| 1056 | WRITE(numout,*) 'Diameter class:', icir |
|---|
| 1057 | WRITE(numout,*) 'Overturning moment multiplier, ', & |
|---|
| 1058 | overturning_moment_multiplier(ipts,ivm) |
|---|
| 1059 | WRITE(numout,*) 'maximum moment for overturning, ', max_overturning_moment(ipts,ivm,icir) |
|---|
| 1060 | WRITE(numout,*) 'maximum moment for stem breakag, ', max_breaking_moment(ipts,ivm,icir) |
|---|
| 1061 | WRITE(numout,*) 'overturning closer,', overturning_moment_closer(ipts,ivm,icir) |
|---|
| 1062 | WRITE(numout,*) 'breaking closer, ', breaking_moment_closer(ipts,ivm,icir) |
|---|
| 1063 | WRITE(numout,*) 'overturning further, ', overturning_moment_further(ipts,ivm,icir) |
|---|
| 1064 | WRITE(numout,*) 'breaking further, ', breaking_moment_further(ipts,ivm,icir) |
|---|
| 1065 | ENDIF |
|---|
| 1066 | !- |
|---|
| 1067 | |
|---|
| 1068 | !! 4. Calculating critical wind speed and damage rate |
|---|
| 1069 | |
|---|
| 1070 | !! 4.1 Calculating critical wind speed in the forest and 10 m above the zero-plane displacement |
|---|
| 1071 | |
|---|
| 1072 | ! +++CHECK+++ |
|---|
| 1073 | ! Ideally ORCHIDEE should use only one surface roughness (z0) and displacement |
|---|
| 1074 | ! height. This is not straightforward because here surface roughness and |
|---|
| 1075 | ! displacement height are calculated as a function of wind speed to account for |
|---|
| 1076 | ! streamlining of the canopy under strong winds. This is done in the subroutine |
|---|
| 1077 | ! streamlining(...) in the calculate_wind_process(). Nevertheless, the subroutine |
|---|
| 1078 | ! streamlining uses the canopy structure from ORCHIDEE and when storm damage |
|---|
| 1079 | ! occurs, the canopy structure is adjusted and this new structure will be used |
|---|
| 1080 | ! in condveg.f90 . The inconsistency is thus between the roughness length |
|---|
| 1081 | ! calculations in condveg.f90 and stomate_windthrow.f90. |
|---|
| 1082 | ! +++++++++++ |
|---|
| 1083 | ! for the edge area |
|---|
| 1084 | ! Critical wind speed for stem breakage in the edge area (m/s) |
|---|
| 1085 | distance_from_force = 0.0 |
|---|
| 1086 | call calculate_wind_process(breaking_moment_closer(ipts,ivm,icir), & |
|---|
| 1087 | current_spacing(icir), canopy_depth(icir), & |
|---|
| 1088 | canopy_breadth(icir), canopy_height(icir), & |
|---|
| 1089 | streamlining_n, streamlining_c, & |
|---|
| 1090 | cws_break_closer(ipts,ivm,icir), & |
|---|
| 1091 | z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir),& |
|---|
| 1092 | distance_from_force) |
|---|
| 1093 | |
|---|
| 1094 | ! Critical wind speed for stem breakage 10 m above |
|---|
| 1095 | ! the zero-plane displacement in the edge area(m/s) |
|---|
| 1096 | call elevate(cws_break_closer(ipts,ivm,icir), & |
|---|
| 1097 | z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir), canopy_height(icir), & |
|---|
| 1098 | cws_break_closer_10(ipts,ivm,icir)) |
|---|
| 1099 | |
|---|
| 1100 | ! Critical wind speed for tree overturning in the edge area (m/s) |
|---|
| 1101 | distance_from_force = 0.0 |
|---|
| 1102 | call calculate_wind_process(overturning_moment_closer(ipts,ivm,icir), & |
|---|
| 1103 | current_spacing(icir), canopy_depth(icir), & |
|---|
| 1104 | canopy_breadth(icir), canopy_height(icir), & |
|---|
| 1105 | streamlining_n, streamlining_c, & |
|---|
| 1106 | cws_overturn_closer(ipts,ivm,icir), & |
|---|
| 1107 | z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir),& |
|---|
| 1108 | distance_from_force) |
|---|
| 1109 | |
|---|
| 1110 | ! Critical wind speed for tree overturning 10 m above |
|---|
| 1111 | ! the zero-plane displacement in the edge area (m/s) |
|---|
| 1112 | call elevate(cws_overturn_closer(ipts,ivm,icir), & |
|---|
| 1113 | z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir), canopy_height(icir), & |
|---|
| 1114 | cws_overturn_closer_10(ipts,ivm,icir)) |
|---|
| 1115 | |
|---|
| 1116 | ! IF the soil at 80 below ground is frozen, the module only allows the stem breakage |
|---|
| 1117 | IF ( soil_temp_daily(ipts) .GT. 273.15 ) THEN |
|---|
| 1118 | ! Calculating damage severity in case any of the critical wind speeds is reached |
|---|
| 1119 | ! in the area around gaps (forest edge area). Whichever of the two damage types |
|---|
| 1120 | ! (overturning, stem breakage) has the lowest critical wind speed, it will be the |
|---|
| 1121 | ! occurring one. |
|---|
| 1122 | IF (cws_overturn_closer_10(ipts,ivm,icir) < cws_break_closer_10(ipts,ivm,icir)) THEN |
|---|
| 1123 | cws_final_closer(ipts,ivm,icir) = cws_overturn_closer_10(ipts,ivm,icir) |
|---|
| 1124 | wind_damage_type_closer(ipts,ivm,icir) = ioverturning |
|---|
| 1125 | ELSE |
|---|
| 1126 | cws_final_closer(ipts,ivm,icir) = cws_break_closer_10(ipts,ivm,icir) |
|---|
| 1127 | wind_damage_type_closer(ipts,ivm,icir) = ibreakage |
|---|
| 1128 | END IF |
|---|
| 1129 | ELSE |
|---|
| 1130 | ! when the soil at 80 cm is frozen, the only damage type is limited to the |
|---|
| 1131 | ! stem_breakage |
|---|
| 1132 | cws_final_closer(ipts,ivm,icir) = cws_break_closer_10(ipts,ivm,icir) |
|---|
| 1133 | wind_damage_type_closer(ipts,ivm,icir) = ibreakage |
|---|
| 1134 | ENDIF |
|---|
| 1135 | |
|---|
| 1136 | ! Based on observations of wind damage in Scotland and how much we had to |
|---|
| 1137 | ! adjust the critical wind speed to get agreement between the wind speed |
|---|
| 1138 | ! above the canopy, the calculated critical wind speed and the presence |
|---|
| 1139 | ! of damage (Hale et al, 2015). I am suggesting the relationship below to |
|---|
| 1140 | ! calculate the percentage of damage at a grid point as a function of the |
|---|
| 1141 | ! meteorological wind speed at 10m above the zero-plane displacement |
|---|
| 1142 | ! (prezhd_wind), and the critical wind speed for damage (crit_wind) |
|---|
| 1143 | ! calculated by ForestGALES. I would set maximum to 80% for the moment |
|---|
| 1144 | ! (sets the maximum total amount of damage = 80%) and the scaling factor |
|---|
| 1145 | ! (s) to 0.8 (this affects the rate at which damage percentage increases |
|---|
| 1146 | ! with increasing wind speed). We can adjust these if the seem to be |
|---|
| 1147 | ! wrong or I find other data that I can use to make a better assessment. |
|---|
| 1148 | ! We apply a sigmodial function to tune the "damage level" by using a |
|---|
| 1149 | ! function which is dependent on the difference between the actual wind |
|---|
| 1150 | ! speed and the final_cws (the lower one of the cws from overturning or |
|---|
| 1151 | ! stem_breakage. |
|---|
| 1152 | wind_damage_rate_closer(ipts,ivm,icir) = max_damage_closer(ivm) * ( & |
|---|
| 1153 | (1./(1.+exp(-((max_wind_speed_storm(ipts)- cws_final_closer(ipts,ivm,icir)) / & |
|---|
| 1154 | sfactor_closer(ivm))))) - & |
|---|
| 1155 | (1./(1.+exp(cws_final_closer(ipts,ivm,icir)/sfactor_closer(ivm)))) ) |
|---|
| 1156 | |
|---|
| 1157 | !-Debug- |
|---|
| 1158 | IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN |
|---|
| 1159 | WRITE(numout,*) 'max_damage_closer:', max_damage_closer(ivm) |
|---|
| 1160 | WRITE(numout,*) 'sfactor_closer:', sfactor_closer(ivm) |
|---|
| 1161 | WRITE(numout,*) 'max_wind_speed_storm (m/s):', max_wind_speed_storm(ipts) |
|---|
| 1162 | WRITE(numout,*) 'cws_closer (m/s)', cws_final_closer(ipts,ivm,icir) |
|---|
| 1163 | WRITE(numout,*) 'wind_damage_rate:', wind_damage_rate_closer(ipts,ivm,icir) |
|---|
| 1164 | ENDIF |
|---|
| 1165 | !- |
|---|
| 1166 | |
|---|
| 1167 | !! 4.2 Calculating critical wind speed in the forest and 10 m above the |
|---|
| 1168 | !! zero-plane displacement for the further forest area |
|---|
| 1169 | |
|---|
| 1170 | ! Whichever of the two damage types (overturning, stem breakage) has the |
|---|
| 1171 | ! lowest critical wind speed, it will be the occurring one. |
|---|
| 1172 | ! Critical wind speed for stem breakage in the inner forest area (m/s) |
|---|
| 1173 | distance_from_force = 0.0 |
|---|
| 1174 | call calculate_wind_process(breaking_moment_further(ipts,ivm,icir), & |
|---|
| 1175 | current_spacing(icir), canopy_depth(icir), & |
|---|
| 1176 | canopy_breadth(icir), canopy_height(icir), & |
|---|
| 1177 | streamlining_n, streamlining_c, & |
|---|
| 1178 | cws_break_further(ipts,ivm,icir), & |
|---|
| 1179 | z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir),& |
|---|
| 1180 | distance_from_force) |
|---|
| 1181 | |
|---|
| 1182 | ! Critical wind speed for stem breakage 10 m above the zero-plane |
|---|
| 1183 | ! displacement in the inner forest area (m/s) |
|---|
| 1184 | call elevate(cws_break_further(ipts,ivm,icir), & |
|---|
| 1185 | z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir), canopy_height(icir), & |
|---|
| 1186 | cws_break_further_10(ipts,ivm,icir)) |
|---|
| 1187 | |
|---|
| 1188 | ! Critical wind speed for tree overturning in the inner forest area (m/s) |
|---|
| 1189 | distance_from_force = 0.0 |
|---|
| 1190 | call calculate_wind_process(overturning_moment_further(ipts,ivm,icir),& |
|---|
| 1191 | current_spacing(icir), canopy_depth(icir), & |
|---|
| 1192 | canopy_breadth(icir), canopy_height(icir), & |
|---|
| 1193 | streamlining_n, streamlining_c, & |
|---|
| 1194 | cws_overturn_further(ipts,ivm,icir),& |
|---|
| 1195 | z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir),& |
|---|
| 1196 | distance_from_force) |
|---|
| 1197 | |
|---|
| 1198 | ! Critical wind speed for tree overturning 10 m above the zero-plane |
|---|
| 1199 | ! displacement in the inner forest area (m/s) |
|---|
| 1200 | call elevate(cws_overturn_further(ipts,ivm,icir), & |
|---|
| 1201 | z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir), canopy_height(icir), & |
|---|
| 1202 | cws_overturn_further_10(ipts,ivm,icir)) |
|---|
| 1203 | |
|---|
| 1204 | ! Whichever of the two damage types (overturning, stem breakage) has the |
|---|
| 1205 | ! lowest critical wind speed, it will be the occurring one. |
|---|
| 1206 | IF (cws_overturn_further_10(ipts,ivm,icir) < cws_break_further_10(ipts,ivm,icir)) THEN |
|---|
| 1207 | cws_final_further(ipts,ivm,icir) = cws_overturn_further_10(ipts,ivm,icir) |
|---|
| 1208 | wind_damage_type_further(ipts,ivm,icir) = ioverturning |
|---|
| 1209 | ELSE |
|---|
| 1210 | cws_final_further(ipts,ivm,icir) = cws_break_further_10(ipts,ivm,icir) |
|---|
| 1211 | wind_damage_type_further(ipts,ivm,icir) = ibreakage |
|---|
| 1212 | END IF |
|---|
| 1213 | |
|---|
| 1214 | ! Wind damage level is considered as a logistic function with sigmoidal |
|---|
| 1215 | ! shape, which is calculated as function of a critical wind speed, an actual |
|---|
| 1216 | ! wind speed, a scaling factor(s)=0.8 and a maximum damage rate(80%)=0.8. |
|---|
| 1217 | wind_damage_rate_further(ipts,ivm,icir) = max_damage_further(ivm) * ( & |
|---|
| 1218 | (1./(1.+exp(-((max_wind_speed_storm(ipts)- cws_final_further(ipts,ivm,icir)) / & |
|---|
| 1219 | sfactor_further(ivm))))) - & |
|---|
| 1220 | (1./(1.+exp(cws_final_further(ipts,ivm,icir)/sfactor_further(ivm)))) ) |
|---|
| 1221 | !-Debug- |
|---|
| 1222 | IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN |
|---|
| 1223 | WRITE(numout,*) 'max_damage_further:', max_damage_further(ivm) |
|---|
| 1224 | WRITE(numout,*) 'sfactor_further:', sfactor_further(ivm) |
|---|
| 1225 | WRITE(numout,*) 'max_wind_speed_storm (m/s):', max_wind_speed_storm(ipts) |
|---|
| 1226 | WRITE(numout,*) 'cws_further (m/s)', cws_final_further(ipts,ivm,icir) |
|---|
| 1227 | WRITE(numout,*) 'wind_damage_rate:', wind_damage_rate_further(ipts,ivm,icir) |
|---|
| 1228 | ENDIF |
|---|
| 1229 | !- |
|---|
| 1230 | |
|---|
| 1231 | !! 5. Determine which trees are killed by wind damage |
|---|
| 1232 | ! We only do the damage for the virtual_circ_class_n GT min_stomate |
|---|
| 1233 | IF (virtual_circ_class_n(ipts,ivm,icir) .GT. min_stomate) THEN |
|---|
| 1234 | |
|---|
| 1235 | !! 5.1. Kill trees in closer area |
|---|
| 1236 | IF (wind_damage_type_closer(ipts,ivm,icir) == ibreakage) THEN |
|---|
| 1237 | kill_break(ipts,ivm,icir) = & |
|---|
| 1238 | wind_damage_rate_closer(ipts,ivm,icir) * circ_class_n(ipts,ivm,icir) * & |
|---|
| 1239 | closer_area_fraction(icir) |
|---|
| 1240 | ! We assume that there is Only one damage type will occur due to |
|---|
| 1241 | ! its lower limiting wind speed |
|---|
| 1242 | kill_uproot(ipts,ivm,icir) = zero |
|---|
| 1243 | ELSE |
|---|
| 1244 | kill_uproot(ipts,ivm,icir) = & |
|---|
| 1245 | wind_damage_rate_closer(ipts,ivm,icir) * & |
|---|
| 1246 | circ_class_n(ipts,ivm,icir) * closer_area_fraction(icir) |
|---|
| 1247 | ! We assume that there is Only one damage type will occur due to |
|---|
| 1248 | ! its lower limited wind speed |
|---|
| 1249 | kill_break(ipts,ivm,icir) = zero |
|---|
| 1250 | ENDIF |
|---|
| 1251 | |
|---|
| 1252 | !! 5.2 Kill trees in further area |
|---|
| 1253 | IF (wind_damage_type_further(ipts,ivm,icir) == ibreakage) THEN |
|---|
| 1254 | ! Total damage = closer + further |
|---|
| 1255 | kill_break(ipts,ivm,icir) = kill_break(ipts,ivm,icir) + & |
|---|
| 1256 | wind_damage_rate_further(ipts,ivm,icir) * & |
|---|
| 1257 | circ_class_n(ipts,ivm,icir) * further_area_fraction(icir) |
|---|
| 1258 | ! We assume that there is only one damage type will occur due to |
|---|
| 1259 | ! its lower limited wind speed |
|---|
| 1260 | kill_uproot(ipts,ivm,icir) = zero |
|---|
| 1261 | ELSE |
|---|
| 1262 | ! Total damage = closer + further |
|---|
| 1263 | kill_uproot(ipts,ivm,icir) = kill_uproot(ipts,ivm,icir) + & |
|---|
| 1264 | wind_damage_rate_further(ipts,ivm,icir) * circ_class_n(ipts,ivm,icir) * & |
|---|
| 1265 | further_area_fraction(icir) |
|---|
| 1266 | ! We assume that there is only one damage type will occur due to |
|---|
| 1267 | ! its lower limited wind speed |
|---|
| 1268 | kill_break(ipts,ivm,icir) = zero |
|---|
| 1269 | ENDIF |
|---|
| 1270 | ELSE |
|---|
| 1271 | |
|---|
| 1272 | ! set no damage for the no vegetation area |
|---|
| 1273 | kill_break(ipts,ivm,icir) = zero |
|---|
| 1274 | kill_uproot(ipts,ivm,icir) = zero |
|---|
| 1275 | END IF |
|---|
| 1276 | |
|---|
| 1277 | ! For the tree height below five meters (7.0m), we assumed that the damage will not |
|---|
| 1278 | ! occur. At the moment, we don't have good solution for the calculation |
|---|
| 1279 | ! of CWS for the small trees, so we just simply keep the small trees. |
|---|
| 1280 | IF ( mean_height(ipts,ivm,icir) .LT. 7.0 ) THEN |
|---|
| 1281 | kill_break(ipts,ivm,icir) = zero |
|---|
| 1282 | kill_uproot(ipts,ivm,icir) = zero |
|---|
| 1283 | ENDIF |
|---|
| 1284 | |
|---|
| 1285 | IF (kill_break(ipts,ivm,icir) .LT. min_stomate) THEN |
|---|
| 1286 | kill_break(ipts,ivm,icir) = zero |
|---|
| 1287 | ENDIF |
|---|
| 1288 | IF (kill_uproot(ipts,ivm,icir) .LT. min_stomate) THEN |
|---|
| 1289 | kill_uproot(ipts,ivm,icir) = zero |
|---|
| 1290 | ENDIF |
|---|
| 1291 | |
|---|
| 1292 | IF(ok_wlsk .AND. ivm .EQ. test_pft .AND. ipts .EQ. test_grid) THEN |
|---|
| 1293 | WRITE(numout,*) 'damage_further',wind_damage_rate_further(ipts,ivm,icir) |
|---|
| 1294 | WRITE(numout,*) 'damage_closer',wind_damage_rate_closer(ipts,ivm,icir) |
|---|
| 1295 | WRITE(numout,*) 'damage_type',wind_damage_type_further(ipts,ivm,icir), & |
|---|
| 1296 | wind_damage_type_closer(ipts,ivm,icir) |
|---|
| 1297 | WRITE(numout,*) 'ccn',circ_class_n(ipts,ivm,icir) |
|---|
| 1298 | ENDIF |
|---|
| 1299 | |
|---|
| 1300 | |
|---|
| 1301 | ! Save mortality to circ_class_cill |
|---|
| 1302 | circ_class_kill(ipts,ivm,icir,forest_managed(ipts,ivm),icut_storm_break) = & |
|---|
| 1303 | kill_break(ipts,ivm,icir) |
|---|
| 1304 | circ_class_kill(ipts,ivm,icir,forest_managed(ipts,ivm),icut_storm_uproot) = & |
|---|
| 1305 | kill_uproot(ipts,ivm,icir) |
|---|
| 1306 | |
|---|
| 1307 | !-Debug- |
|---|
| 1308 | IF (printlev_loc .GT. 4 .AND. ipts==test_grid .AND. ivm==test_pft) THEN |
|---|
| 1309 | WRITE(numout,*) 'End of Section 4: Calculating critical wind speed and damage rate' |
|---|
| 1310 | WRITE(numout,*) 'critical wind speed further, ', cws_final_further(ipts,ivm,icir) |
|---|
| 1311 | WRITE(numout,*) 'critical wind speed closer, ', cws_final_closer(ipts,ivm,icir) |
|---|
| 1312 | WRITE(numout,*) 'wind damage type further, ', wind_damage_type_further(ipts,ivm,icir) |
|---|
| 1313 | WRITE(numout,*) 'wind damage type closer, ', wind_damage_type_closer(ipts,ivm,icir) |
|---|
| 1314 | WRITE(numout,*) 'forest managed type, ', forest_managed(ipts,ivm) |
|---|
| 1315 | WRITE(numout,*) 'circ_class_n(ipts,ivm,icir),', circ_class_n(ipts,ivm,icir) |
|---|
| 1316 | WRITE(numout,*) 'circ_class_kill(ipts,ivm,icir,:,icut_storm_break,', & |
|---|
| 1317 | circ_class_kill(ipts,ivm,icir,forest_managed(ipts,ivm),icut_storm_break) |
|---|
| 1318 | WRITE(numout,*) 'circ_class_kill(ipts,ivm,icir,:,icut_storm_uproot,', & |
|---|
| 1319 | circ_class_kill(ipts,ivm,icir,forest_managed(ipts,ivm),icut_storm_uproot) |
|---|
| 1320 | WRITE(numout,*) 'circ_class_kill, for icut_1:9_:' |
|---|
| 1321 | WRITE(numout,*) circ_class_kill(ipts,ivm,icir,forest_managed(ipts,ivm),1:9) |
|---|
| 1322 | ENDIF |
|---|
| 1323 | !- |
|---|
| 1324 | |
|---|
| 1325 | ! Copy some values for the output files |
|---|
| 1326 | canopy_depth_out(ipts,ivm,icir) = canopy_depth(icir) |
|---|
| 1327 | tree_h_edge_out(ipts,ivm,icir) = tree_heights_from_edge(icir) |
|---|
| 1328 | mean_gap_f_out(ipts,ivm,icir) = mean_gap_factor(icir) |
|---|
| 1329 | gap_size_out(ipts,ivm) = gap_size |
|---|
| 1330 | |
|---|
| 1331 | END DO ! loop ncirc |
|---|
| 1332 | |
|---|
| 1333 | END DO ! loop nvm |
|---|
| 1334 | |
|---|
| 1335 | END DO ! loop npts |
|---|
| 1336 | |
|---|
| 1337 | !! 5. Writing out history variables |
|---|
| 1338 | |
|---|
| 1339 | ! Write the CWS to stomate_history_file |
|---|
| 1340 | DO icir = 1,ncirc |
|---|
| 1341 | !CWS at further area for each diameter class |
|---|
| 1342 | WRITE(var_name,'(A,I3.3)') 'CWS_FURTHER_',icir |
|---|
| 1343 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1344 | cws_final_further(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1345 | ENDDO |
|---|
| 1346 | ! |
|---|
| 1347 | DO icir =1,ncirc |
|---|
| 1348 | !CWS at closer area for each diameter class |
|---|
| 1349 | WRITE(var_name,'(A,I3.3)') 'CWS_CLOSER_',icir |
|---|
| 1350 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1351 | cws_final_closer(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1352 | ENDDO |
|---|
| 1353 | ! |
|---|
| 1354 | DO icir =1,ncirc |
|---|
| 1355 | !CWS at closer area in each diameter class for stem breakage |
|---|
| 1356 | WRITE(var_name,'(A,I3.3)') 'CWS_CLOSER_BK_',icir |
|---|
| 1357 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1358 | cws_break_closer_10(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1359 | ENDDO |
|---|
| 1360 | ! |
|---|
| 1361 | DO icir =1,ncirc |
|---|
| 1362 | !CWS at closer area in each diameter class for overturning |
|---|
| 1363 | WRITE(var_name,'(A,I3.3)') 'CWS_CLOSER_OV_',icir |
|---|
| 1364 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1365 | cws_overturn_closer_10(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1366 | ENDDO |
|---|
| 1367 | ! |
|---|
| 1368 | DO icir =1,ncirc |
|---|
| 1369 | !CWS at further area in each diameter class for stem break |
|---|
| 1370 | WRITE(var_name,'(A,I3.3)') 'CWS_FURTHER_BK_',icir |
|---|
| 1371 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1372 | cws_break_further_10(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1373 | ENDDO |
|---|
| 1374 | ! |
|---|
| 1375 | DO icir =1,ncirc |
|---|
| 1376 | !CWS at further area in each diameter class for overturning |
|---|
| 1377 | WRITE(var_name,'(A,I3.3)') 'CWS_FURTHER_OV_',icir |
|---|
| 1378 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1379 | cws_overturn_further_10(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1380 | ENDDO |
|---|
| 1381 | ! |
|---|
| 1382 | DO icir =1,ncirc |
|---|
| 1383 | !CWS at closer area in each diameter class for stem breakage |
|---|
| 1384 | WRITE(var_name,'(A,I3.3)') 'CWS_TOP_CLO_BK_',icir |
|---|
| 1385 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1386 | cws_break_closer(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1387 | ENDDO |
|---|
| 1388 | ! |
|---|
| 1389 | DO icir =1,ncirc |
|---|
| 1390 | !CWS at closer area in each diameter class for overturning |
|---|
| 1391 | WRITE(var_name,'(A,I3.3)') 'CWS_TOP_CLO_OV_',icir |
|---|
| 1392 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1393 | cws_overturn_closer(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1394 | ENDDO |
|---|
| 1395 | ! |
|---|
| 1396 | DO icir =1,ncirc |
|---|
| 1397 | !CWS at further area in each diameter class for stem break |
|---|
| 1398 | WRITE(var_name,'(A,I3.3)') 'CWS_TOP_FUR_BK_',icir |
|---|
| 1399 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1400 | cws_break_further(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1401 | ENDDO |
|---|
| 1402 | ! |
|---|
| 1403 | DO icir =1,ncirc |
|---|
| 1404 | !CWS at further area in each diameter class for overturning |
|---|
| 1405 | WRITE(var_name,'(A,I3.3)') 'CWS_TOP_FUR_OV_',icir |
|---|
| 1406 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1407 | cws_overturn_further(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1408 | ENDDO |
|---|
| 1409 | ! |
|---|
| 1410 | DO icir =1,ncirc |
|---|
| 1411 | !D_H_ratio for each diameter class |
|---|
| 1412 | WRITE(var_name,'(A,I3.3)') 'D_H_RATIO_',icir |
|---|
| 1413 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1414 | d_h_ratio(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1415 | ENDDO |
|---|
| 1416 | ! |
|---|
| 1417 | DO icir =1,ncirc |
|---|
| 1418 | ! MEAN_HEIGHT for each diameter class |
|---|
| 1419 | WRITE(var_name,'(A,I3.3)') 'MEAN_HEIGHT_',icir |
|---|
| 1420 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1421 | mean_height(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1422 | ENDDO |
|---|
| 1423 | ! |
|---|
| 1424 | DO icir =1,ncirc |
|---|
| 1425 | ! MEAN_DBH for each diameter class |
|---|
| 1426 | WRITE(var_name,'(A,I3.3)') 'MEAN_DBH_',icir |
|---|
| 1427 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1428 | mean_dbh(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1429 | ENDDO |
|---|
| 1430 | ! |
|---|
| 1431 | DO icir =1,ncirc |
|---|
| 1432 | ! V_IND virtual individual tree numbers for each diameter class |
|---|
| 1433 | WRITE(var_name,'(A,I3.3)') 'V_IND_',icir |
|---|
| 1434 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1435 | virtual_circ_class_n(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1436 | ENDDO |
|---|
| 1437 | ! |
|---|
| 1438 | DO icir =1,ncirc |
|---|
| 1439 | ! G_CLOSER gustiness for closer area in each diameter class |
|---|
| 1440 | WRITE(var_name,'(A,I3.3)') 'G_CLOSER_',icir |
|---|
| 1441 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1442 | gust_factor_closer(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1443 | ENDDO |
|---|
| 1444 | ! |
|---|
| 1445 | DO icir =1,ncirc |
|---|
| 1446 | ! G_FURTHER gustiness for further area in each diameter class |
|---|
| 1447 | WRITE(var_name,'(A,I3.3)') 'G_FURTHER_',icir |
|---|
| 1448 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1449 | gust_factor_further(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1450 | ENDDO |
|---|
| 1451 | ! |
|---|
| 1452 | DO icir =1,ncirc |
|---|
| 1453 | ! G_FURTHER gustiness for further area in each diameter class |
|---|
| 1454 | WRITE(var_name,'(A,I3.3)') 'STEM_MASS_',icir |
|---|
| 1455 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1456 | stem_mass(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1457 | ENDDO |
|---|
| 1458 | ! |
|---|
| 1459 | DO icir =1,ncirc |
|---|
| 1460 | ! roughness length |
|---|
| 1461 | WRITE(var_name,'(A,I3.3)') 'WIND_Z0_',icir |
|---|
| 1462 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1463 | z0_wind(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1464 | ENDDO |
|---|
| 1465 | ! |
|---|
| 1466 | DO icir =1,ncirc |
|---|
| 1467 | ! displacement height |
|---|
| 1468 | WRITE(var_name,'(A,I3.3)') 'WIND_D_',icir |
|---|
| 1469 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1470 | zhd_wind(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1471 | ENDDO |
|---|
| 1472 | ! |
|---|
| 1473 | DO icir =1,ncirc |
|---|
| 1474 | ! canopy depth |
|---|
| 1475 | WRITE(var_name,'(A,I3.3)') 'CANOPY_DEPTH_',icir |
|---|
| 1476 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1477 | canopy_depth_out(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1478 | ENDDO |
|---|
| 1479 | ! |
|---|
| 1480 | DO icir =1,ncirc |
|---|
| 1481 | ! tree heights from edge |
|---|
| 1482 | WRITE(var_name,'(A,I3.3)') 'TREE_H_EDGE_',icir |
|---|
| 1483 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1484 | tree_h_edge_out(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1485 | ENDDO |
|---|
| 1486 | ! |
|---|
| 1487 | DO icir =1,ncirc |
|---|
| 1488 | ! mean gap factor |
|---|
| 1489 | WRITE(var_name,'(A,I3.3)') 'MEAN_GAP_F_',icir |
|---|
| 1490 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1491 | mean_gap_f_out(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1492 | ENDDO |
|---|
| 1493 | ! |
|---|
| 1494 | ! gap size |
|---|
| 1495 | WRITE(var_name,'(A)') 'GAP_SIZE' |
|---|
| 1496 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1497 | gap_size_out(:,:), npts*nvm, horipft_index) |
|---|
| 1498 | ! |
|---|
| 1499 | DO icir =1,ncirc |
|---|
| 1500 | ! edge factor |
|---|
| 1501 | WRITE(var_name,'(A,I3.3)') 'EDGE_FACTOR_',icir |
|---|
| 1502 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1503 | edge_factor_out(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1504 | ENDDO |
|---|
| 1505 | ! |
|---|
| 1506 | DO icir =1,ncirc |
|---|
| 1507 | ! maximum moment for overturning |
|---|
| 1508 | WRITE(var_name,'(A,I3.3)') 'MAX_OV_MOMENT_',icir |
|---|
| 1509 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1510 | max_overturning_moment(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1511 | ENDDO |
|---|
| 1512 | ! |
|---|
| 1513 | DO icir =1,ncirc |
|---|
| 1514 | ! maximum moment breakage |
|---|
| 1515 | WRITE(var_name,'(A,I3.3)') 'MAX_BK_MOMENT_',icir |
|---|
| 1516 | CALL histwrite_p(hist_id_stomate, var_name, itime, & |
|---|
| 1517 | max_breaking_moment(:,:,icir), npts*nvm, horipft_index) |
|---|
| 1518 | ENDDO |
|---|
| 1519 | ! |
|---|
| 1520 | ! Daily Maximum wind speed |
|---|
| 1521 | CALL histwrite_p(hist_id_stomate, 'DAILY_MAX_WIND', itime, & |
|---|
| 1522 | wind_speed_daily, npts, hori_index) |
|---|
| 1523 | |
|---|
| 1524 | CALL xios_orchidee_send_field("CWS_FURTHER", cws_final_further) |
|---|
| 1525 | CALL xios_orchidee_send_field("CWS_CLOSWER", cws_final_closer) |
|---|
| 1526 | CALL xios_orchidee_send_field("CWS_CLOSER_BK", cws_break_closer_10) |
|---|
| 1527 | CALL xios_orchidee_send_field("CWS_CLOSER_OV", cws_overturn_closer_10) |
|---|
| 1528 | CALL xios_orchidee_send_field("CWS_FURTHER_BK", cws_break_further_10) |
|---|
| 1529 | CALL xios_orchidee_send_field("CWS_FURTHER_OV", cws_overturn_further_10) |
|---|
| 1530 | CALL xios_orchidee_send_field("GAP_AREA",gap_area_save(:,:,1)) |
|---|
| 1531 | |
|---|
| 1532 | CALL xios_orchidee_send_field("DAILY_MAX_WIND", wind_speed_daily) |
|---|
| 1533 | CALL xios_orchidee_send_field("V_IND", virtual_circ_class_n) |
|---|
| 1534 | CALL xios_orchidee_send_field("KILL_STORM_BK", kill_break) |
|---|
| 1535 | CALL xios_orchidee_send_field("KILL_STORM_OV", kill_uproot) |
|---|
| 1536 | |
|---|
| 1537 | END SUBROUTINE wind_damage |
|---|
| 1538 | |
|---|
| 1539 | |
|---|
| 1540 | !! ============================================================================================= |
|---|
| 1541 | !! SUBROUTINE : calculate_wind_process |
|---|
| 1542 | !! |
|---|
| 1543 | !>\BRIEF : This subroutine calculates the breaking and overturning critical wind |
|---|
| 1544 | !! speeds |
|---|
| 1545 | !! |
|---|
| 1546 | !! DESCRIPTION : None |
|---|
| 1547 | !! |
|---|
| 1548 | !! RECENT CHANGE(S): None |
|---|
| 1549 | !! |
|---|
| 1550 | !! RETURN VALUE : |
|---|
| 1551 | !! |
|---|
| 1552 | !! REFERENCES : Gardiner et al. 2010, hale et al. 2012, Raupach 1994. |
|---|
| 1553 | !! |
|---|
| 1554 | !! FLOWCHART : None |
|---|
| 1555 | !! |
|---|
| 1556 | !_ ============================================================================================= |
|---|
| 1557 | |
|---|
| 1558 | SUBROUTINE calculate_wind_process(critical_moment,& |
|---|
| 1559 | current_spacing, canopy_depth, canopy_breadth, canopy_height, & |
|---|
| 1560 | streamlining_n, streamlining_c, & |
|---|
| 1561 | guess_wind_speed, z0_wind, zhd_wind, distance_from_force) |
|---|
| 1562 | |
|---|
| 1563 | IMPLICIT NONE |
|---|
| 1564 | |
|---|
| 1565 | !! 0. Variable and parameter declaration |
|---|
| 1566 | |
|---|
| 1567 | !! 0.1 Input variables |
|---|
| 1568 | REAL, INTENT(in) :: current_spacing |
|---|
| 1569 | REAL, INTENT(in) :: canopy_depth |
|---|
| 1570 | REAL, INTENT(in) :: canopy_breadth |
|---|
| 1571 | REAL, INTENT(in) :: canopy_height |
|---|
| 1572 | REAL, INTENT(in) :: critical_moment |
|---|
| 1573 | REAL, INTENT(in) :: streamlining_n |
|---|
| 1574 | REAL, INTENT(in) :: streamlining_c |
|---|
| 1575 | REAL, INTENT(in) :: distance_from_force |
|---|
| 1576 | |
|---|
| 1577 | !! 0.2 Output variables |
|---|
| 1578 | |
|---|
| 1579 | REAL, INTENT(out) :: guess_wind_speed !! A guess for the critical wind speed |
|---|
| 1580 | REAL, INTENT(out) :: z0_wind |
|---|
| 1581 | REAL, INTENT(out) :: zhd_wind |
|---|
| 1582 | !! 0.3 Modified variables |
|---|
| 1583 | |
|---|
| 1584 | !! 0.4 Local variables |
|---|
| 1585 | INTEGER :: loop_count !! Counting the iteration |
|---|
| 1586 | INTEGER :: max_count=20 !! Max iteration times set to 20 iterations |
|---|
| 1587 | REAL :: delta !! Difference between the guessed |
|---|
| 1588 | !! and the calculated wind speed |
|---|
| 1589 | REAL :: wind_precision !! How close we should |
|---|
| 1590 | !! get with the iterations |
|---|
| 1591 | REAL :: lambda !! |
|---|
| 1592 | REAL :: lambda_capital |
|---|
| 1593 | REAL :: gamma_solved |
|---|
| 1594 | REAL :: psih |
|---|
| 1595 | REAL :: critical_moment_cws |
|---|
| 1596 | |
|---|
| 1597 | REAL, PARAMETER :: air_density = 1.2250 |
|---|
| 1598 | REAL, PARAMETER :: f_crown_mass = 1.136 |
|---|
| 1599 | REAL, PARAMETER :: f_knots = 1.000 |
|---|
| 1600 | |
|---|
| 1601 | !! 0.5 Variables used in calculate_force subroutine and calculate_bending_moment subroutine |
|---|
| 1602 | !! (Both subroutines are not used in the current code) |
|---|
| 1603 | REAL :: gammaSolved |
|---|
| 1604 | REAL :: d |
|---|
| 1605 | REAL :: force_on_tree |
|---|
| 1606 | REAL :: height_of_force |
|---|
| 1607 | REAL :: bending_moment |
|---|
| 1608 | REAL :: porosity |
|---|
| 1609 | REAL :: crit_wind_speed |
|---|
| 1610 | |
|---|
| 1611 | !_ ============================================================================================= |
|---|
| 1612 | |
|---|
| 1613 | !! 1. Subroutine |
|---|
| 1614 | |
|---|
| 1615 | ! This bit belongs to the original form of doing the iterations. For more info, see the notes below. |
|---|
| 1616 | ! guess_wind_speed = 64.0 |
|---|
| 1617 | ! delta = guess_wind_speed / 2.0 |
|---|
| 1618 | |
|---|
| 1619 | guess_wind_speed = 25.0 |
|---|
| 1620 | delta = guess_wind_speed / 2.0 |
|---|
| 1621 | wind_precision = 0.01 |
|---|
| 1622 | loop_count = 0 |
|---|
| 1623 | DO WHILE (delta > wind_precision) |
|---|
| 1624 | CALL streamlining(guess_wind_speed, current_spacing, & |
|---|
| 1625 | canopy_breadth, canopy_depth, canopy_height, & |
|---|
| 1626 | streamlining_n, streamlining_c, gamma_solved, zhd_wind, z0_wind) |
|---|
| 1627 | |
|---|
| 1628 | ! Calculate critical wind speed to compare with wind speeds used for calculating |
|---|
| 1629 | ! streamlining in the iterative calculation of the final critical wind speed) |
|---|
| 1630 | crit_wind_speed = (gamma_solved / current_spacing) * (critical_moment / & |
|---|
| 1631 | (f_crown_mass * air_density * (zhd_wind - distance_from_force) ))**0.5 |
|---|
| 1632 | |
|---|
| 1633 | ! We use the absolute error as the criteria for stopping the iteration |
|---|
| 1634 | ! loop |
|---|
| 1635 | delta = abs( guess_wind_speed - crit_wind_speed ) |
|---|
| 1636 | |
|---|
| 1637 | guess_wind_speed = crit_wind_speed |
|---|
| 1638 | |
|---|
| 1639 | loop_count = loop_count +1 |
|---|
| 1640 | ! Add the maximum iteration times condition for CWS calculation |
|---|
| 1641 | IF (loop_count .GE. max_count) THEN |
|---|
| 1642 | WRITE(numout,*) 'WARANNING from stomate_windthrow.f90: The maximum iteration times for CWS claculation!!' |
|---|
| 1643 | EXIT |
|---|
| 1644 | ENDIF |
|---|
| 1645 | END DO |
|---|
| 1646 | |
|---|
| 1647 | END SUBROUTINE calculate_wind_process |
|---|
| 1648 | |
|---|
| 1649 | |
|---|
| 1650 | !! =========================================================================================== |
|---|
| 1651 | !! SUBROUTINE : streamlining |
|---|
| 1652 | !! |
|---|
| 1653 | !>\BRIEF : This subroutine calculates the streamlining of tree crowns. Streamlining is |
|---|
| 1654 | !! the change of shape of the crowns due to wind. |
|---|
| 1655 | !! |
|---|
| 1656 | !! DESCRIPTION : None |
|---|
| 1657 | !! |
|---|
| 1658 | !! RECENT CHANGE(S): None |
|---|
| 1659 | !! |
|---|
| 1660 | !! RETURN VALUE : Surface roughness and zero-plane displacement among other variables that |
|---|
| 1661 | !! are used in the other subroutines |
|---|
| 1662 | !! REFERENCES : |
|---|
| 1663 | !! Raupach (1994) Simplified expressions for vegetation roughness length and |
|---|
| 1664 | !! zero-plane displacement as functions of canopy height and area index, Boundary-Layer Meteorology |
|---|
| 1665 | !! October 1994, Volume 71, Issue 1, pp 211-216, doi: 10.1007/BF00709229 |
|---|
| 1666 | !! |
|---|
| 1667 | !! FLOWCHART : None |
|---|
| 1668 | !! |
|---|
| 1669 | !_ ============================================================================================ |
|---|
| 1670 | |
|---|
| 1671 | SUBROUTINE streamlining(guess_wind_speed, current_spacing, & |
|---|
| 1672 | canopy_breadth, canopy_depth, a_mean_height, & |
|---|
| 1673 | streamlining_n, streamlining_c, gamma_solved, zhd_wind, z0_wind) |
|---|
| 1674 | |
|---|
| 1675 | IMPLICIT NONE |
|---|
| 1676 | |
|---|
| 1677 | !! 0. Variable and parameter declaration |
|---|
| 1678 | |
|---|
| 1679 | !! 0.1 Input variables |
|---|
| 1680 | |
|---|
| 1681 | REAL, INTENT(in) :: guess_wind_speed !! the guess wind speed (m s{-1}) |
|---|
| 1682 | REAL, INTENT(in) :: current_spacing !! the average spacing between trees (m) |
|---|
| 1683 | REAL, INTENT(in) :: canopy_breadth !! the maximum width of the canopy (m) |
|---|
| 1684 | REAL, INTENT(in) :: canopy_depth !! the length of the live tree crown (m) |
|---|
| 1685 | REAL, INTENT(in) :: a_mean_height !! the mean canopy height (m) |
|---|
| 1686 | REAL, INTENT(in) :: streamlining_n !! streamlining parameter n (-) see Eq. A12 in Hale et al 2015 |
|---|
| 1687 | REAL, INTENT(in) :: streamlining_c !! streamlining parameter c (-) see Eq. A12 in Hale et al 2015 |
|---|
| 1688 | |
|---|
| 1689 | !! 0.2 Output variables |
|---|
| 1690 | |
|---|
| 1691 | REAL, INTENT(out) :: gamma_solved !! a function for solving surface roughness (-) |
|---|
| 1692 | REAL, INTENT(out) :: zhd_wind !! zero plane of displacement height (m) |
|---|
| 1693 | REAL, INTENT(out) :: z0_wind !! surface roughness (m) |
|---|
| 1694 | |
|---|
| 1695 | !! 0.3 Modified variables |
|---|
| 1696 | |
|---|
| 1697 | !! 0.4 Local variables |
|---|
| 1698 | |
|---|
| 1699 | REAL :: lambda !! roughness density, the frontal |
|---|
| 1700 | !! area of roughness elements |
|---|
| 1701 | REAL :: lambda_capital !! canopy area index |
|---|
| 1702 | REAL :: psih !! stability correction function for heat transfer (-) |
|---|
| 1703 | REAL :: porosity !! the effective drag ?? (-) |
|---|
| 1704 | !! see equation A12 in the supplementary material of |
|---|
| 1705 | !! Hale et al (2015) |
|---|
| 1706 | |
|---|
| 1707 | !! The aerodynamic parameter values used here are based on the reference of Raupach (1994) |
|---|
| 1708 | !! --- TEMP ---- |
|---|
| 1709 | !! Need to be discussed for the consistency in the ORCHIDEE |
|---|
| 1710 | !! ------------ |
|---|
| 1711 | REAL, PARAMETER :: c_displacement = 7.5 |
|---|
| 1712 | REAL, PARAMETER :: c_surface = 0.003 |
|---|
| 1713 | REAL, PARAMETER :: c_drag = 0.3 |
|---|
| 1714 | REAL, PARAMETER :: c_roughness = 2.0 |
|---|
| 1715 | REAL, PARAMETER :: lambda_capital_max = 0.6 |
|---|
| 1716 | REAL, PARAMETER :: wind_streamlining_max = 25.0 |
|---|
| 1717 | REAL, PARAMETER :: wind_streamlining_min = 10.0 |
|---|
| 1718 | !_ ============================================================================================= |
|---|
| 1719 | |
|---|
| 1720 | !! 1. Subroutine |
|---|
| 1721 | |
|---|
| 1722 | IF (guess_wind_speed > wind_streamlining_max) THEN |
|---|
| 1723 | ! Above 25 m/s, tree crowns cannot streamline any further. |
|---|
| 1724 | ! The porosity set to its minimum value 2.35*25^-0.51 ~ 0.46 |
|---|
| 1725 | porosity = streamlining_c * wind_streamlining_max**(-streamlining_n) |
|---|
| 1726 | ELSE IF ( guess_wind_speed < wind_streamlining_min ) THEN |
|---|
| 1727 | ! Below 10 m/s, tree crowns do not streamline and porosity set to its |
|---|
| 1728 | ! maximum value 2.35*10^-0.51 ~ 0.73 |
|---|
| 1729 | porosity = streamlining_c * wind_streamlining_min**(-streamlining_n) |
|---|
| 1730 | ELSE |
|---|
| 1731 | ! Wind speed between 10 m/s and 25 m/s, tree crown do the streamlining |
|---|
| 1732 | porosity = streamlining_c * guess_wind_speed**(-streamlining_n) |
|---|
| 1733 | END IF |
|---|
| 1734 | |
|---|
| 1735 | ! lambda = (canopy_breadth/2.0) * canopy_depth * porosity / current_spacing **2.0 |
|---|
| 1736 | ! Here, we directly calculated the mean width of canopy from the value of "canopy_breadth" |
|---|
| 1737 | ! The shape of tree crown is not rhomboid anymore. In the ORCHIDEE-CAN |
|---|
| 1738 | ! the shape of tree crown is circle, so the frontal area should be |
|---|
| 1739 | ! calculated as (pi/4 * canopy_breadth * canopy_depth). |
|---|
| 1740 | lambda = canopy_breadth * canopy_depth * porosity / current_spacing **2.0 |
|---|
| 1741 | |
|---|
| 1742 | lambda_capital = lambda |
|---|
| 1743 | |
|---|
| 1744 | IF (lambda_capital > lambda_capital_max) THEN |
|---|
| 1745 | ! Gamma is needed for the calculation of the roughness length. |
|---|
| 1746 | gamma_solved = 1.0 / ((c_surface + c_drag * lambda_capital_max / 2.0)**0.5) |
|---|
| 1747 | ELSE |
|---|
| 1748 | gamma_solved = 1.0 / ((c_surface + c_drag * lambda_capital / 2.0)**0.5) |
|---|
| 1749 | END IF |
|---|
| 1750 | |
|---|
| 1751 | ! Calculation of the zero-plane displacement: |
|---|
| 1752 | ! see equations from A7 to A11 in the supplementary material of Hale et al. (2015) |
|---|
| 1753 | zhd_wind = a_mean_height * (1.0 - ((1.0 - EXP(-(c_displacement * & |
|---|
| 1754 | lambda_capital)**0.5)) / (c_displacement * lambda_capital)**0.5)) |
|---|
| 1755 | |
|---|
| 1756 | ! Calculation of the aerodynamic roughness length: |
|---|
| 1757 | ! psih = "Profile influence function" at h (height of the roughness elements). |
|---|
| 1758 | psih = LOG(c_roughness) -1.0 + c_roughness**(-1.0) |
|---|
| 1759 | |
|---|
| 1760 | z0_wind = (a_mean_height - zhd_wind) * EXP(psih - ct_karman * gamma_solved) |
|---|
| 1761 | |
|---|
| 1762 | END SUBROUTINE streamlining |
|---|
| 1763 | |
|---|
| 1764 | |
|---|
| 1765 | !! |
|---|
| 1766 | !============================================================================================= |
|---|
| 1767 | !! SUBROUTINE : elevate |
|---|
| 1768 | !! |
|---|
| 1769 | !>\BRIEF : This subroutine converts the critical wind speeds |
|---|
| 1770 | !calculated above to the |
|---|
| 1771 | !! wind speed 10 meters above the zero-plane displacement |
|---|
| 1772 | !! |
|---|
| 1773 | !! DESCRIPTION : None |
|---|
| 1774 | !! |
|---|
| 1775 | !! RECENT CHANGE(S): None |
|---|
| 1776 | !! |
|---|
| 1777 | !! RETURN VALUE : |
|---|
| 1778 | !! |
|---|
| 1779 | !! REFERENCES : |
|---|
| 1780 | !! |
|---|
| 1781 | !! FLOWCHART : None |
|---|
| 1782 | !! |
|---|
| 1783 | !_ |
|---|
| 1784 | !============================================================================================= |
|---|
| 1785 | |
|---|
| 1786 | SUBROUTINE elevate(uh_speed, z0_wind, zhd_wind, mean_height,elevate_result) |
|---|
| 1787 | |
|---|
| 1788 | IMPLICIT NONE |
|---|
| 1789 | |
|---|
| 1790 | !! 0. Variable and parameter declaration |
|---|
| 1791 | |
|---|
| 1792 | !! 0.1 Input variables |
|---|
| 1793 | |
|---|
| 1794 | REAL, INTENT(in) :: uh_speed |
|---|
| 1795 | REAL, INTENT(in) :: z0_wind |
|---|
| 1796 | REAL, INTENT(in) :: zhd_wind |
|---|
| 1797 | REAL, INTENT(in) :: mean_height |
|---|
| 1798 | |
|---|
| 1799 | !! 0.2 Output variables |
|---|
| 1800 | |
|---|
| 1801 | REAL, INTENT(out) :: elevate_result |
|---|
| 1802 | |
|---|
| 1803 | !! 0.3 Modified variables |
|---|
| 1804 | |
|---|
| 1805 | !! 0.4 Local variables |
|---|
| 1806 | |
|---|
| 1807 | !_ |
|---|
| 1808 | !============================================================================================== |
|---|
| 1809 | |
|---|
| 1810 | !! 1. Subroutine |
|---|
| 1811 | elevate_result = uh_speed * log(10.0 / z0_wind) / & |
|---|
| 1812 | log((mean_height - zhd_wind) / z0_wind) |
|---|
| 1813 | |
|---|
| 1814 | END SUBROUTINE elevate |
|---|
| 1815 | |
|---|
| 1816 | |
|---|
| 1817 | !============================================================================================== |
|---|
| 1818 | ! Below subroutines are not used for now |
|---|
| 1819 | !============================================================================================== |
|---|
| 1820 | |
|---|
| 1821 | !! ============================================================================================ |
|---|
| 1822 | !! SUBROUTINE : calculate_force |
|---|
| 1823 | !! |
|---|
| 1824 | !>\BRIEF : This subroutine calculates the force affecting the tree. |
|---|
| 1825 | !! |
|---|
| 1826 | !! DESCRIPTION : None |
|---|
| 1827 | !! |
|---|
| 1828 | !! RECENT CHANGE(S): None |
|---|
| 1829 | !! |
|---|
| 1830 | !! RETURN VALUE : |
|---|
| 1831 | !! |
|---|
| 1832 | !! REFERENCES : |
|---|
| 1833 | !! |
|---|
| 1834 | !! FLOWCHART : None |
|---|
| 1835 | !! |
|---|
| 1836 | !_ ============================================================================================ |
|---|
| 1837 | |
|---|
| 1838 | SUBROUTINE calculate_force(guess_wind_speed, a_gamma_solved, & |
|---|
| 1839 | a_zhd_wind, current_spacing, a_mean_height, & |
|---|
| 1840 | force_on_tree, height_of_force) |
|---|
| 1841 | |
|---|
| 1842 | IMPLICIT NONE |
|---|
| 1843 | |
|---|
| 1844 | !! 0. Variable and parameter declaration |
|---|
| 1845 | |
|---|
| 1846 | !! 0.1 Input variables |
|---|
| 1847 | |
|---|
| 1848 | REAL, INTENT(in) :: guess_wind_speed |
|---|
| 1849 | REAL, INTENT(in) :: a_gamma_solved |
|---|
| 1850 | REAL, INTENT(in) :: a_zhd_wind |
|---|
| 1851 | REAL, INTENT(in) :: current_spacing |
|---|
| 1852 | REAL, INTENT(in) :: a_mean_height |
|---|
| 1853 | !! 0.2 Output variables |
|---|
| 1854 | |
|---|
| 1855 | REAL, INTENT(out) :: force_on_tree |
|---|
| 1856 | REAL, INTENT(out) :: height_of_force |
|---|
| 1857 | |
|---|
| 1858 | !! 0.3 Modified variables |
|---|
| 1859 | |
|---|
| 1860 | !! 0.4 Local variables |
|---|
| 1861 | REAL,PARAMETER :: air_density = 1.2226 !!This should call from constant.f90 (kg m^{-3}) |
|---|
| 1862 | !_ ============================================================================================= |
|---|
| 1863 | |
|---|
| 1864 | !! 1. Subroutine |
|---|
| 1865 | |
|---|
| 1866 | force_on_tree = air_density * (guess_wind_speed * & |
|---|
| 1867 | current_spacing / a_gamma_solved)**2.0 |
|---|
| 1868 | |
|---|
| 1869 | !! --- TEMP --- |
|---|
| 1870 | !! I checked with the reference the force of wind can be calculated as |
|---|
| 1871 | !! air_density*ustar^{2}*D^{2} as the gamma is defined as U/ustar , D is the current spacing |
|---|
| 1872 | !! the force of wind should be as |
|---|
| 1873 | !! force_on_tree = air_density * (U/gamma)^{2} * D^{2} |
|---|
| 1874 | !! ------------ |
|---|
| 1875 | height_of_force = a_zhd_wind |
|---|
| 1876 | |
|---|
| 1877 | ! The maximum of the zero-plane displacement height is locked at 0.8 tree height. |
|---|
| 1878 | ! IF (height_of_force > 0.8 * a_mean_height) THEN |
|---|
| 1879 | ! height_of_force = 0.8 * a_mean_height |
|---|
| 1880 | ! END IF |
|---|
| 1881 | |
|---|
| 1882 | END subroutine calculate_force |
|---|
| 1883 | |
|---|
| 1884 | |
|---|
| 1885 | !! ============================================================================================= |
|---|
| 1886 | !! SUBROUTINE : calculate_bending_moment |
|---|
| 1887 | !! |
|---|
| 1888 | !>\BRIEF : Subroutine to calculate the bending moment |
|---|
| 1889 | !! (to compare with critical moments in the iterative calculation of the critical |
|---|
| 1890 | !! wind speeds) |
|---|
| 1891 | !! |
|---|
| 1892 | !! DESCRIPTION : None |
|---|
| 1893 | !! |
|---|
| 1894 | !! RECENT CHANGE(S): None |
|---|
| 1895 | !! |
|---|
| 1896 | !! RETURN VALUE : |
|---|
| 1897 | !! |
|---|
| 1898 | !! REFERENCES : |
|---|
| 1899 | !! |
|---|
| 1900 | !! FLOWCHART : None |
|---|
| 1901 | !! |
|---|
| 1902 | !_ ============================================================================================= |
|---|
| 1903 | |
|---|
| 1904 | SUBROUTINE calculate_bending_moment (a_force_on_tree, a_height_of_force, bending_moment) |
|---|
| 1905 | |
|---|
| 1906 | IMPLICIT NONE |
|---|
| 1907 | |
|---|
| 1908 | !! 0. Variable and parameter declaration |
|---|
| 1909 | |
|---|
| 1910 | !! 0.1 Input variables |
|---|
| 1911 | |
|---|
| 1912 | REAL, INTENT(in) :: a_force_on_tree |
|---|
| 1913 | REAL, INTENT(in) :: a_height_of_force |
|---|
| 1914 | |
|---|
| 1915 | !! 0.2 Output variables |
|---|
| 1916 | |
|---|
| 1917 | REAL, INTENT(out) :: bending_moment |
|---|
| 1918 | |
|---|
| 1919 | !! 0.3 Modified variables |
|---|
| 1920 | |
|---|
| 1921 | !! 0.4 Local variables |
|---|
| 1922 | REAL :: f_crown_mass |
|---|
| 1923 | |
|---|
| 1924 | !_ ============================================================================================= |
|---|
| 1925 | |
|---|
| 1926 | !! 1. Subroutine |
|---|
| 1927 | |
|---|
| 1928 | bending_moment = a_force_on_tree * a_height_of_force * f_crown_mass |
|---|
| 1929 | |
|---|
| 1930 | END SUBROUTINE calculate_bending_moment |
|---|
| 1931 | |
|---|
| 1932 | |
|---|
| 1933 | END MODULE stomate_windthrow |
|---|
