! ================================================================================================================================= ! MODULE : stomate_wind ! ! CONTACT : yiyingchen@gate.sinica.edu.tw ! ! LICENCE : IPSL (2006). This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC ! !>\BRIEF : Calculates critical wind speed for tree damage for each half-hourly period of the simulation. ! This wind speed can be compared to the actual wind speed to decide whether a tree damage occurs. !! !!\streamlining_n DESCRIPTION: None !! !! RECENT CHANGE(S): None !! !! SVN : !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/stomate_windthrow.f90 $ !! $Date: 2015 -01-22 13:41:13 +0100 (mer., 22 janv. 2014) $ !! $Revision: 1831 $ !! \streamlining_n !_ ================================================================================================================================ MODULE stomate_windthrow ! modules used: USE ioipsl_para USE grid USE xios_orchidee USE stomate_data USE constantes USE constantes_soil USE pft_parameters USE pft_parameters_var USE function_library, ONLY: wood_to_dia, Nmax, wood_to_qmdia,wood_to_height, & wood_to_cv, cc_to_lai, cc_to_biomass, biomass_to_lai IMPLICIT NONE ! private & public routines PRIVATE PUBLIC wind_damage !!$ wind_clear, streamlining, calculate_force, calculate_bending_moment, & !!$ calculate_wind_process, elevate LOGICAL, SAVE :: firstcall_windthrow = .TRUE. !! First call !$OMP THREADPRIVATE(firstcall_windthrow) INTEGER(i_std), SAVE :: printlev_loc !! Local level of text output for current module !$OMP THREADPRIVATE(printlev_loc) CONTAINS !! ================================================================================================================================ !! SUBROUTINE : wind_clear !! !>\BRIEF Set the firstcall flag to .TRUE. and activate initialization !! !_ ================================================================================================================================ SUBROUTINE wind_clear firstcall_windthrow = .TRUE. END SUBROUTINE wind_clear !! ================================================================================================================================ !! SUBROUTINE : wind_damage !! !>\BRIEF This subroutine calculates the critical wind speeds for both uprooting and stem breakage !! for each half-hourly period of the simulation. The algorithm follows Barry Gardiner's wind !! damage risk model called GALES (Hale et al. 2015). !! !! DESCRIPTION : The purpose of this module is to describe the actual sensitivity of a given PFT to the wind that is present !! at the same pixel. For each circumference class in each PFT, two critical wind speed values (CWS) are calculated in every !! half-hourly time step by this subroutine, following the approach of differentiating two types of wind damage of trees: !! - uprooting (the whole tree is removed from the soil together with its root plate) !! - stem breakage (the root plate and the bottom part of the stem remains attached to the ground). !! Whichever threshold (i.e. CWS) is reached first for a tree, the according damage type occurs in the model. !! !! The module completes the following tasks: !! 1. Calculates... !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): :: Critical wind speed for stem breakage 10 m above the zero-plane displacement (m/s): U_Break_10 !! Critical wind speed for tree overturn 10 m above the zero-plane displacement (m/s): U_Overturn_10 !! !! REFERENCES : !! - Hale, S., Gardiner, B., Nicoll, B., Taylor, P., and Pizzirani, S. (2015). Comparison and validation of three versions !! of a forest wind risk model. Environmental Modelling and Software. In Press. !! ... !! !! FLOWCHART : None !! \streamlining_n !_ ================================================================================================================================ SUBROUTINE wind_damage (npts, nlevels_tot, circ_class_biomass, & veget_max, circ_class_n, wind_speed_daily, plant_status, & circ_class_kill, gap_area_save, forest_managed, & soil_temp_daily, root_profile, max_wind_speed_storm, count_storm,& is_storm) IMPLICIT NONE !! 0. Variable declarations !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: npts !! Domain size @tex $(unitless)$ @endtex INTEGER(i_std), INTENT(in) :: nlevels_tot !! Total level numbers in photosynthesis INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: forest_managed !! forest management flag (is the forest !! being managed?) REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass in of an individual tree in each circ class !! @tex $(m^{-2})$ @endtex REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of trees within each circumference REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT on the ground !! May sum to less than unity if the pixel has nobio !! area. (unitless, 0-1) REAL(r_std), DIMENSION(:,:), INTENT(in) :: plant_status !! Growth and phenological status of the plant REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: gap_area_save !! Gap area created by more than 30% of the basal area loss !! in the last 5 years !! Dimension(npts,nvm,wind_years) REAL(r_std), DIMENSION(:), INTENT(in) :: soil_temp_daily !! Daily maximum soil temperature at 0.8 meter below ground (K) REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: root_profile !! Normalized root mass/length fraction in each soil layer !! (0-1, unitless) !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std), DIMENSION(:,:,:,:,:), & INTENT(inout) :: circ_class_kill !! Number of trees within a circ that needs !! to be killed @tex $(ind m^{-2})$ @endtex REAL(r_std), DIMENSION(:), INTENT(inout) :: wind_speed_daily !! Daily maximum wind speed at 2 meter (ms-1) REAL(r_std), DIMENSION(:), INTENT(inout) :: max_wind_speed_storm !! Daily maximum wind speed at 2 meter (ms-1) during a storm event INTEGER(i_std), DIMENSION(:), INTENT(inout) :: count_storm !! number of day after a storm LOGICAL, DIMENSION(:), INTENT(inout) :: is_storm !! are we in a storm event ? !! 0.4 Local variables CHARACTER(30) :: var_name !! To store variable names for I/O INTEGER(i_std) :: islm,ipts,ivm !! Index for the loops INTEGER(i_std) :: icir,i,iele, iyear !! Index for the loops INTEGER(i_std), DIMENSION(npts) :: soil_type !! Soil types: !! 1. Free-draining mineral soils; !! 2. Gleyed mineral soils; !! 3. Peaty mineral soils; and !! 4. Deep Peats REAL(r_std), DIMENSION(npts,nvm,ncirc) :: virtual_circ_class_n !! Virtual number of trees within each circumference REAL(r_std), DIMENSION(ncirc,60) :: section_height !! The height of the bottom of the 1-m sections !! of the mean tree @tex $(m)$ @endtex. 60 is the !! max number of sections thus max tree height is 60m REAL(r_std), DIMENSION(ncirc,60) :: branch_width !! Width of canopy) of the 1-m sections of the mean !! tree @tex $(m)$ @endtex 60 is the max number of !! sections thus max tree height is 60m REAL(r_std), DIMENSION(npts,nvm,ncirc) :: mean_dbh !! Diameter at breast height (dbh) of the mean tree !! of the circumference-class @tex $(m)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: mean_height !! Height of the mean tree of the given dbh-class !! @tex $(m)$ @endtex REAL(r_std), DIMENSION(ncirc) :: lai !! Leaf area index @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: canopy_breadth !! Canopy breadth (width of canopy) @tex $(m)$ @endtex REAL(r_std), DIMENSION(ncirc) :: canopy_depth !! Canopy depth (height of canopy) @tex $(m)$ @endtex REAL(r_std), DIMENSION(ncirc) :: canopy_height !! The height of the starting point of the canopy from !! below @tex $(m)$ @endtex REAL(r_std), DIMENSION(ncirc) :: mid_canopy !! The middle height of the canopy @tex $(m)$ @endtex REAL(r_std), DIMENSION(ncirc) :: mean_dbh_height !! The height of the diameter (by default it is at !! breast height = 1.3 m) @tex $(m)$ @endtex REAL(r_std), DIMENSION(ncirc) :: snow_mass !! Weight of the snow on the tree @tex $(kg)$ @endtex REAL(r_std), DIMENSION(ncirc) :: diameter_c2 !! Variable for calculating diameters !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: diameter_canopy_base !! Diameter at the base of the canopy !! @tex $(m)$ @endtex REAL(r_std), DIMENSION(ncirc) :: diameter_c1 !! Variable for calculating diameters !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: crown_depth !! Crown depth of the tree in each circumference class REAL(r_std), DIMENSION(ncirc) :: crown_width !! Crown width of the mean tree (m) REAL(r_std), DIMENSION(ncirc) :: crown_volume !! Crown volume of the mean tree @tex $(m^{3})$ @endtex REAL(r_std), DIMENSION(ncirc) :: crown_mass !! Mass of the mean tree crown @tex $(kg)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: stem_mass !! Mass of the mean tree stem @tex $(kg)$ @endtex REAL(r_std), DIMENSION(ncirc) :: current_spacing !! Spacing between trees in the given dbh-class !! @tex $(m)$ @endtexi INTEGER(i_std), DIMENSION(npts,nvm) :: rooting_depth !! Rooting depth of the mean tree (three categories; !! 1: shallow, 2: deep, 3:average !! See constantes_soil_var.f90 for the definitions REAL(r_std) :: streamlining_c !! Streamlining parameter. @tex $(unitless)$ @endtex. !! Streamlining is the change of shape of the crowns !! due to wind. REAL(r_std) :: streamlining_n !! Streamlining parameter. @tex $(unitless)$ @endtex. !! Streamlining is the change of shape of the crowns !! due to wind. REAL(r_std) :: area_timber_removals_5_years !! Area of the total timber removal in the !! previous 5 year @tex $(m^{2})$ @endtex REAL(r_std), DIMENSION(ncirc) :: tree_heights_from_edge !! Width of the forest edge area, i.e. area_closer (in number of tree heights) !! @tex $(unitless)$ @endtex REAL(r_std) :: n_gap !! The number of gaps with an area of clear_cut_max @tex $(unitless)$ @endtex REAL(r_std) :: length_side_gap !! The length of one side of the square-shaped gap @tex $(m)$ @endtex REAL(r_std), DIMENSION(ncirc) :: area_around_gap !! The framing edge area around the gap with different gustiness !! @tex $(m^{2})$ @endtex REAL(r_std), DIMENSION(ncirc) :: area_total_closer !! The total area of forest edge areas affected by wind @tex $(m^{2})$ @endtex REAL(r_std), DIMENSION(ncirc) :: area_total_further !! The total area of forest without edge areas @tex $(m^{2})$ @endtex REAL(r_std), DIMENSION(ncirc) :: ratio_spacing_height !! Ratio of tree height and tree spacing @tex $(unitless)$ @endtex REAL(r_std) :: gap_size !! The length (width of the cross-section) of the gap @tex $(m)$ @endtex REAL(r_std), DIMENSION(ncirc) :: mean_gap_factor !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: max_gap_factor !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: term_a !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: term_b !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: old_gust_closer !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: old_gust_further !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: old_gust_edge !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: old_gust_edge_gap !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: new_gust_closer !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: new_gust_further !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: new_gust_edge !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: new_gust_edge_gap_closer !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: new_gust_edge_gap_further !! Variable used for calculating gustiness !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(ncirc) :: mean_bm_gust_factor_closer !! Mean bending moment of the forest edge area !! (for calculating the edge factor) !! @tex $(Nm/kg)$ @endtex REAL(r_std), DIMENSION(ncirc) :: mean_bm_gust_factor_further !! Mean bending moment of the inner forest area (for calculating !! the edge factor) @tex $(Nm/kg)$ @endtex REAL(r_std), DIMENSION(ncirc) :: edge_factor !! The ratio of the mean bending moment in the forest edge area !! and mean bending moment in the inner forest area !! @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: gust_factor_closer !! Gust factor for calculating the critical wind speed of the forest !! edge area @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: gust_factor_further !! Gust factor for calculating the critical wind speed of the inner !! forest area @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: g_closer !! Gustiness of the forest edge area @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: g_further !! Gustiness of the inner forest area @tex $(unitless)$ @endtex REAL(r_std), DIMENSION(npts,nvm) :: overturning_moment_multiplier!! Overturning moment multiplier representing the effects of species, !! soil type, soil depths and senescence @tex $(Nm/kg)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: max_overturning_moment !! The maximum overturning moment the mean tree can withstand !! @tex $(Nm/kg)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: max_breaking_moment !! The maximum breaking moment the mean tree can withstand !! @tex $(Nm/kg)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: overturning_moment_closer !! The maximum overturning moment the mean tree can withstand in the !! forest edge area with gustiness taken into account @tex $(Nm/kg)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: overturning_moment_further !! The maximum overturning moment the mean tree can withstand in the !! inner forest area with gustiness taken into account @tex $(Nm/kg)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: breaking_moment_closer !! The maximum breaking moment the mean tree can withstand in the !! forest edge area with gustiness taken into account @tex $(Nm/kg)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: breaking_moment_further !! The maximum breaking moment the mean tree can withstand in the !! inner forest area with gustiness taken into account @tex $(Nm/kg)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_break_closer !! Critical wind speed for stem breakage in the forest edge area !! @tex $(m/s; mean hourly wind speed)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_break_closer_10 !! Critical wind speed for stem breakage in the forest edge area !! @tex $(m/s; mean hourly wind speed)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_overturn_closer !! Critical wind speed for overturning in the forest edge area !! @tex $(m/s; mean hourly wind speed)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_overturn_closer_10 !! Critical wind speed for overturning in the forest edge area !! @tex $(m/s; mean hourly wind speed)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_final_closer !! The critical wind speed of the final form of damage (breakage / uprooting) !! in the forest edge area @tex $(m/s; mean hourly wind speed)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: d_h_ratio !! Array for storing the stand spacing and tree height ratio for output (-) REAL(r_std), DIMENSION(npts) :: wind_speed_actual !! The actual wind speed over the pixel !! @tex $(m/s; mean hourly wind speed)$ @endtex LOGICAL :: wind_damage_type_uprooting_closer !! The final form of damage (breakage / uprooting) in the forest edge area REAL(r_std), DIMENSION(npts,nvm,ncirc) :: wind_damage_rate_closer !! Damage rate (%) in the forest edge area in case of damage !! @tex $(0-1)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_break_further !! Critical wind speed for stem breakage in the inner forest area !! @tex $(m/s; mean hourly wind speed)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_break_further_10 !! Critical wind speed for stem breakage in the inner forest area !! @tex $(m/s; mean hourly wind speed)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_overturn_further !! Critical wind speed for overturning in the inner forest area !! @tex $(m/s; mean hourly wind speed)$ @endtex !Jina: Seems it if for 10 m above the z0 but has the same description. REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_overturn_further_10 !! Critical wind speed for overturning in the inner forest area !! @tex $(m/s; mean hourly wind speed)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: cws_final_further !! The critical wind speed of the final form of damage (breakage / uprooting) !! in the inner forest area @tex $(m/s; mean hourly wind speed)$ @endtex LOGICAL :: wind_damage_type_uprooting_further !! The final form of damage (breakage / uprooting) in the inner forest area REAL(r_std), DIMENSION(npts,nvm,ncirc) :: wind_damage_rate_further !! Damage rate (%) in the forest edge area in case !! of damage in the inner forest area @tex $(0-1)$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: z0_wind !! surface roughness calculated from the subroutine streamlining() REAL(r_std), DIMENSION(npts,nvm,ncirc) :: zhd_wind !! zero plane displacement height calculated from the subroutine streamlining(m) INTEGER(i_std), DIMENSION(npts,nvm,ncirc) :: wind_damage_type_closer !! Integer store the wind damage type for closer area ibreakage=1, ioverturning=2 INTEGER(i_std), DIMENSION(npts,nvm,ncirc) :: wind_damage_type_further !! Integer store the wind damage type for further area REAL(r_std), DIMENSION(npts,nvm,ncirc) :: kill_break !! Biomass loss through stem breakage from storm !! @tex $(ind m^{-2})$ @endtex REAL(r_std), DIMENSION(npts,nvm,ncirc) :: kill_uproot !! Number of trees killed through uproot from storm !! @tex $(ind m^{-2})$ @endtex REAL(r_std) :: crown_top !! Temporal variable for storing the tallest tree crown height in different PFTs REAL(r_std) :: crown_bottom !! Temporal variable for storing the smallest tree crown height in different PFTs REAL(r_std) :: distance_from_force !! The distance from the force was set to 1.3 for breakage and 0.0 for overturning REAL(r_std), DIMENSION(ncirc) :: closer_area_fraction !! Closer area fraction REAL(r_std), DIMENSION(ncirc) :: further_area_fraction !! Further area fraction REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: biomass ! Temp output Variables for GALES model intercomparision REAL(r_std), DIMENSION(npts,nvm,ncirc) :: canopy_depth_out REAL(r_std), DIMENSION(npts,nvm,ncirc) :: tree_h_edge_out REAL(r_std), DIMENSION(npts,nvm,ncirc) :: mean_gap_f_out REAL(r_std), DIMENSION(npts,nvm) :: gap_size_out REAL(r_std), DIMENSION(npts,nvm,ncirc) :: edge_factor_out REAL(r_std) :: roots_below !! temporary variable containing the share of roots below 80 cm depth (0-1, unitless) !_ ============================================================================================================================================================================ IF (firstcall_windthrow) THEN !! Initialize local printlev printlev_loc=get_printlev('stomate_windthrow') ! The code is specific for CRU-NCEP. If a different forcing ! is used, the Max Wind Ratio parameters should be adjusted CALL ipslerr_p(2,'Wind throw module assumes that CRU-NCEP forcing is used',& '','','') END IF IF (printlev_loc>=2) WRITE(numout,*) 'Entering stomate_windthrow.f90' !! 0. Initialized the values of variables wind_damage_rate_closer(:,:,:) = zero wind_damage_rate_further(:,:,:) = zero cws_final_further(:,:,:) = zero cws_final_closer(:,:,:) = zero d_h_ratio(:,:,:) = zero kill_uproot = zero kill_break = zero cws_break_further_10 = zero cws_break_closer_10 = zero cws_overturn_further_10 = zero cws_overturn_closer_10 = zero virtual_circ_class_n = zero !-Debug- IF (printlev_loc .GT. 4) THEN DO ipts = 1, npts WRITE(numout,*) 'Daily maximum windspeed passed from stomate to stomate_winthrow:', wind_speed_daily(ipts) WRITE(numout,*) 'Daily maximum soil temperature at 80cm below ground:', soil_temp_daily(ipts) ENDDO ENDIF !- !! 1. Getting the input variables in the required form for the actual subroutine ! calculate the total biomass biomass(:,:,:,:) = cc_to_biomass(npts,nvm,circ_class_biomass,circ_class_n) !! 1.1 Convert 6-hr CRU-NCEP into half-hourly wind speed ! The max_wind_ratio parameters were determined by relating fluxnet ! half-hourly data to 6h CRU-NCEP data. For this reason the ! parameters should be considered specific to the CRU-NCEP forcing. ! The max_wind_ratio parameters were also found to depend on the ! spatial aggregation, hence, they are specific to the 0.5 degree ! CRU-NCEP data (see Chen et al for details). ! The following relationship to convert CRU-NCEP into half-hourly ! wind speeds was found: !wind_speed_daily(:) = -5.922*wind_speed_daily(:) & ! + 2.051*wind_speed_daily(:)**2. & ! - 0.191*wind_speed_daily(:)**3. & ! + 0.006*wind_speed_daily(:)**4. !WHERE (wind_speed_daily(:) .LE. zero) ! wind_speed_daily(:) = zero !ENDWHERE ! In case this convertion did not result in the observed wind ! damage, we kept a linear tuning parameter. Ideally, no tuning ! should be applied and therefore daily_max_tune should be 1.0 !wind_speed_daily(:) = wind_speed_daily(:) * daily_max_tune ! Storms do not last for exactly 24 hours. ! Storms that pass overnight would be accounted for twice ! unless we add criteria to tell ORCHIDEE this is just a ! single storm. The value wind_speed_damage_thr ! (set to 20.0 m/s) corresponds to the wind speed ! threshold above which is_storm flag is set the TRUE ! and nb_days_storm (set to 5.0 days) correspond to the ! number of days at which the max wind speed is less than the ! threshold which means that the storm is over and hence ! is_storm flag to FALSE. As a consequence the windstorm ! damage is calculated nb_days_storm (set to 5 days) ! after the end of a storm. One main consequences of this ! new approach is that damages due to wind is only accounted ! when the flag is_storm is set to TRUE. WHERE (wind_speed_daily(:) .GT. & wind_speed_storm_thr .AND. .NOT. is_storm(:)) is_storm(:) = .TRUE. count_storm(:) = 0 ENDWHERE DO ipts = 1, npts if(is_storm(ipts)) then max_wind_speed_storm(ipts)= & max(max_wind_speed_storm(ipts), & wind_speed_daily(ipts)) if(wind_speed_daily(ipts) .LT. & wind_speed_storm_thr) then count_storm(ipts) = count_storm(ipts) + 1 endif if(count_storm(ipts) .LT. nb_days_storm) then cycle endif count_storm(ipts) = 0 is_storm(:) = .FALSE. else cycle endif !! 1.3 Defining the soil type ! This sub-routine uses 4 types of soil for calculating the ! critical wind speed, so 12 different soil types present in ! other modules of ORCHIDEE need to be grouped into these ! 4 generic types. ! +++CHECK+++ ! For the moment only 3 soil types are distinguished. ! Needs to be checked with all 12 USDA soil types. ! If more soil type are used it will become possible ! to make use of the different soils distinguished in ! the windthrow module. The code could then look as ! follow: !!$ SELECT CASE() !!$ CASE (x:y) !!$ soil_type(ipts) = ifreedrainage !!$ CASE (x:y) !!$ soil_type(ipts) = igleyed !!$ CASE (x:y) !!$ soil_type(ipts) = ipeaty !!$ CASE (x:y) !!$ soil_type(ipts) = ipeat !!$ CASE DEFAULT !!$ CALL ipslerr_p(3, 'stomate_windthrow.f90',& !!$ 'soil type is not defined','','') !!$ END SELECT ! For the moment all soils are assumed to be freely draining soil_type(ipts) = ifree_draining ! +++++++++++++ DO ivm = 2, nvm ! Initilaize area_timber_removals_5_years = zero n_gap = zero ! Don't do the calculations if there is no vegetation ! or the vegetation is bare soil, grass or crop IF ((veget_max(ipts,ivm).LT.min_stomate) .OR. .NOT.is_tree(ivm)) THEN CYCLE ENDIF !! 1.4 Calculating the rooting depth ! If root density (proportion of root mass) is higher ! than 40% below 80 cm below the soil surface, then we ! treat rooting depth as "deep". Otherwise, it's "shallow". ! Note that when using the static root profile deep vs shallow ! will be the same for all pixels belonhing to the same PFT ! beacuse it solely depends on humcste. The threshold of 40% ! was chosen to ensure that all trees are shallow rooting with ! the static root profile. IF (grnd_80.EQ.nslm) THEN ! Handle the case where the deepest soil layer is less ! than 80 cm. All roots must then be above 80 cm. roots_below = un ELSE ! Calculate the share of roots below 80 cm roots_below = zero DO islm = grnd_80,nslm roots_below = roots_below + root_profile(ipts,ivm,islm,istruc) END DO END IF IF (roots_below > 0.40) THEN rooting_depth(ipts,ivm) = ideep ELSE rooting_depth(ipts,ivm) = ishallow END IF !--Debug IF (printlev_loc >= 4) THEN WRITE(numout,*) 'rooting depth ', rooting_depth(ipts,ivm) ENDIF !! 1.5 Set crown parameters ! To set crown parameters leaf mass needs to be available IF (SUM(circ_class_biomass(ipts,ivm,:,ileaf,icarbon)) .GT. min_stomate) THEN ! There is a canopy streamlining_c = streamlining_c_leaf(ivm) streamlining_n = streamlining_n_leaf(ivm) ELSE ! If the tree is leafless, the following parameters ! become different. The weight of leaves will be absent. ! Streamlining changes as the surface of leaves is absent. streamlining_c = streamlining_c_leafless(ivm) streamlining_n = streamlining_n_leafless(ivm) END IF !-Debug- IF(printlev_loc .GT. 4 .AND. ivm .EQ. test_pft)THEN WRITE(numout,*) 'Initialisation of critical windspeed' WRITE(numout,*) 'ipts, ivm, ',ipts, ivm WRITE(numout,*) 'soil_type(ipts), ',soil_type(ipts) WRITE(numout,*) 'rooting_depth_type,(1:deep, 2:shallow,3:average), ',& rooting_depth(ipts,ivm) WRITE(numout,*) 'streamlining_c, ', streamlining_c WRITE(numout,*) 'streamlining_n, ', streamlining_n ENDIF !- !! 1.6 Calculate stand characteristics ! Calculating diameter at breast height (dbh) in m of the ! mean tree of all circ classes. mean_dbh(ipts,ivm,:) = wood_to_dia(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm) ! Calculating the height in m of the mean tree of the given ! dbh-class using mean_height(ipts,ivm,:) = wood_to_height(& circ_class_biomass(ipts,ivm,:,:,icarbon),ivm) ! Crown_volume is calculated as the woody biomass of an individual tree crown_volume(:) = wood_to_cv( circ_class_biomass(ipts,ivm,:,:,icarbon), & circ_class_n(ipts,ivm,:),ivm) ! Crown_width is calculated as(crown_vloume * 6.0/pi)^(1/3) crown_width(:) = ((crown_volume(:)*6.0)/pi)**(1.0/3.0) ! ORCHIDEE assumes that the crowns are spherical. Hence, width ! breath and depth are all the same. Set canopy_breadth as the ! crown_width for all circumference classes canopy_breadth(:) = crown_width(:) canopy_depth(:) = crown_width(:) canopy_height(:) = mean_height(ipts,ivm,:) ! Calculate gaps and the surrounding areas ! IMPORTANT: AS WE DO NOT HAVE INFORMATION ON STAND EDGES IN THE ! ORCHIDEE SIMULATIONS, A POSSIBLE PROXY TO USE HERE IS A VARIABLE ! THAT REPRESENTS EARLIER ! GAP AREAS IN THE PREVIOUS 5 YEARS. AFTER ! SUCH A PERIOD, THE FOREST EDGE IS ASSUMED TO GET ACCLIMATIZED ! TO THE NEW CIRCUMSTANCES, AND IT'S VULNERABILITY DIMINISHES ! (see Persson, 1975; Valinger and Fridman, 2011) ! Note that gap area for the current year is not calculated yet (which is ! calculated at the end of the year in stomate_season), so the gap ! created t-1 to t-5 is applied. ! Calculating the total area of timber removals in the previous ! five years (set by the variable ::wind_years area_timber_removals_5_years = zero DO iyear = 1, wind_years area_timber_removals_5_years = area_timber_removals_5_years + & gap_area_save(ipts,ivm,iyear) END DO ! Number of gaps with the area of the maximum allowed clearfelling: n_gap = area_timber_removals_5_years / clear_cut_max ! The length of the sides of the squared-shaped gaps: length_side_gap = SQRT(clear_cut_max) !-Debug- IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN WRITE(numout,*) 'crown_volumn, ', crown_volume(:) WRITE(numout,*) 'crown_width, ', crown_width(:) WRITE(numout,*) 'crown_depth, ', canopy_depth(:) WRITE(numout,*) 'canopy_height, ', mean_height(:,:,:) ENDIF !- ! loop over the circ class DO icir = 1, ncirc !-Debug- IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN WRITE(numout,*) 'biomass, isapabove, ', biomass(ipts,ivm,isapabove,icarbon) WRITE(numout,*) 'biomass, iheartabove, ', biomass(ipts,ivm,iheartabove,icarbon) WRITE(numout,*) 'circ_class_biomass, isapabove, ', & circ_class_biomass(ipts,ivm,icir,isapabove,icarbon) WRITE(numout,*) 'circ_class_biomass, iheatabove, ', & circ_class_biomass(ipts,ivm,icir,iheartabove,icarbon) ENDIF !- !! 1.6 Calculating the current tree spacing in the PFT ! The current tree spacing between tree stems are calculated ! from the number of stems per hectare. circ_class_n is the ! number of stems per m2; it's converted to stems per hectare ! with m2_to_ha which has a value of 10000. ! The virtual tree spacing assumes that the LAI remains ! constant and the calculates how many trees of a given ! circumference class would be needed to result in the set LAI. ! Because lai changes throughout the season because of leaf ! fall, wood mass was used instead. Because LAI is related to ! wood mass (see allocation). The principle remains the same IF ( (circ_class_biomass(ipts,ivm,icir,isapabove,icarbon) + & circ_class_biomass(ipts,ivm,icir,iheartabove,icarbon)) >= zero ) THEN virtual_circ_class_n(ipts,ivm,icir) = & (biomass(ipts,ivm,isapabove,icarbon) + & biomass(ipts,ivm,iheartabove,icarbon)) / & (circ_class_biomass(ipts,ivm,icir,isapabove,icarbon) + & circ_class_biomass(ipts,ivm,icir,iheartabove,icarbon)) ELSE virtual_circ_class_n(ipts,ivm,icir) = zero ENDIF ! Now the spacing is calculated based on the stem density ! based on the virtual circ class to beeter deal with ! unmanaged forest and heterogenous canopy. IF (virtual_circ_class_n(ipts,ivm,icir) .GT. min_stomate) THEN ! The original equation is current_spacing(icir) = 100.0 / & ! SQRT(m2_to_ha * virtual_circ_class_n(ipts,ivm,icir)) ! given that m2_to_ha = 10000 it reduces to: current_spacing(icir) = 1.0 / SQRT(virtual_circ_class_n(ipts,ivm,icir)) ELSE ! In principle, we should never enter this condition. The ! spacing between 2 trees equals the size of the pixel. It is ! thus empty. We defined this case such that if we do end here, ! a very large cws will be calculated and the code should not ! suffer from divide by zero. current_spacing(icir) = SQRT(area(ipts)) WRITE(numout,*) 'WARNING from WINDFALL module:' WRITE(numout,*) 'the current spacing set to pixel spacing:', current_spacing(icir) CALL ipslerr_p(2,'Wind throw module',& 'current spacing was set to pixel dimensions',& 'it is therefore assumed to be empty','') ENDIF !! 1.7 Calculate green stem mass ! Stem mass is calculated as the aboveground woody biomass ! minus the branches. circ_class_biomass is the biomass of ! an individual tree (gC/tree) ! Here, we convert the unit to (kgC/tree) by dividing kilo_to_unit (1000) stem_mass(ipts,ivm,icir) = (circ_class_biomass(ipts,ivm,icir,isapabove,icarbon) + & circ_class_biomass(ipts,ivm,icir,iheartabove,icarbon)) * & (un - branch_ratio(ivm)) / kilo_to_unit !-Debug- IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN WRITE(numout,*) 'stem_mass,',stem_mass(ipts,ivm,icir) WRITE(numout,*) 'green_density,', green_density(ivm) WRITE(numout,*) 'pipe_sensity,', pipe_density(ivm) WRITE(numout,*) 'kilo to unit,', kilo_to_unit WRITE(numout,*) 'virtual_circ_class', virtual_circ_class_n(ipts,ivm,icir) WRITE(numout,*) 'sla,',sla(ivm) END IF !- ! Wind throw occurs on living trees so the stem mass should be ! converted from dry biomass in gC to wet biomass in kg to be used ! in the GALES equations. Stem_mass thus converted as follows: ! (green_density / (pipe_density/1000)) from dry wood mass ! (kgC/tree) to wet wood mass (Kg/trees) stem_mass(ipts,ivm,icir) = stem_mass(ipts, ivm,icir) * & green_density(ivm)/(pipe_density(ivm)/kilo_to_unit) !! 1.8 Calculating leaf area index ! We need to use the virtual biomass in each circumference class ! to calculate lai, and tree_heights_from_edge was calculated as ! function of lai. Use the function biomass_to_lai to account for ! a dynamic lai calculation. lai(icir) = biomass_to_lai(circ_class_biomass(ipts,ivm,icir,ileaf,icarbon)*& virtual_circ_class_n(ipts,ivm,icir),ivm) !-Debug- IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN WRITE(numout,*) 'check spacing and circ_number in a diameter class' WRITE(numout,*) 'virtual_tree_numbers:', & virtual_circ_class_n(ipts,ivm,icir), ' for circ class:', icir WRITE(numout,*) 'current spacing: ',& current_spacing(icir),'mean diameter: ', mean_dbh(ipts,ivm,icir) WRITE(numout,*) 'area:', area(ipts) WRITE(numout,*) 'ipts, ',ipts,'ivm, ',ivm,'icir, ',icir WRITE(numout,*) 'mean dbh, ', mean_dbh(ipts,ivm,icir) WRITE(numout,*) 'mean height, ', mean_height(ipts,ivm,icir) WRITE(numout,*) 'stem_mass, ', stem_mass(ipts,ivm,icir) WRITE(numout,*) 'crown_volume, ',crown_volume(icir) WRITE(numout,*) 'crown_mass, ',crown_mass(icir) WRITE(numout,*) 'crown_width, ',crown_width(icir) WRITE(numout,*) 'lai, ',lai(icir) WRITE(numout,*) 'virtual_circ_n, ', virtual_circ_class_n(ipts,ivm,icir) WRITE(numout,*) 'real_circ_n, ', circ_class_n(ipts,ivm,icir) WRITE(numout,*) 'virtual_spacing, ',current_spacing(icir) WRITE(numout,*) 'real_spacing, ', SQRT(1.0/circ_class_n(ipts,ivm,icir)) WRITE(numout,*) 'End of Section 1: Calculated canopy structure & & in difffernt circumferences for each PFT' ENDIF ! !! 2. Calculating the gustiness of wind in the PFTs !! 2.1 Calculate gaps and the surrounding areas ! IMPORTANT: AS WE DO NOT HAVE INFORMATION ON STAND EDGES IN THE ! ORCHIDEE SIMULATIONS, A POSSIBLE PROXY TO USE HERE IS A VARIABLE ! THAT REPRESENTS EARLIER HARVESTS IN THE PREVIOUS 5 YEARS. AFTER ! SUCH A PERIOD, THE FOREST EDGE IS ASSUMED TO GET ACCLIMATIZED ! TO THE NEW CIRCUMSTANCES, AND IT'S VULNERABILITY DIMINISHES ! (see Persson, 1975; Valinger and Fridman, 2011) ! There is a big increase in gustiness after a certain distance from ! the edge of the forest, hence we divide the area of the PFT into two ! areas: ! - the area which is further than X tree heights from any edge that ! was created in the previous five years (area_total_further) ! - the area which is closer than X tree heights from any edge that ! was created in the previous five years (area_total_closer). ! Of course, these smaller areas depend on the fragmentation of the ! gaps (created in the last five years). As we have no information ! on this, we assume a fixed gap size representing the most likely ! size of clear cuts (clear_cut_max) and that the ! gaps are square-shaped. Then, we divide the area of harvest in the ! previous five years (area_timber_removals_5_years) by clear_cut_max. ! Hence, from the number of gaps we can calculate the area_total_closer. ! As each circumference class will have its forest area divided, we ! calculate two gustiness values and hence, two critical wind speeds. ! The value of X depends on the leaf area index of the PFT and is ! calculated as follows: ! tree_heights_from_edge is defined as "X/h", X sets to 28/LAI and h ! is set to meant_tree_height. ! This comes from the equation of calculating canopy penetration depth ! in Pan et al. 2017?. Ying Pan work on the large eddy simulation and ! is now at Pennsilvania State University ! Here, lai(icir) for each diameter class is calculated by virtual_class_n ! a threshold "0.1" is used to avoid the issue of divided by small lai < 0.1 ! for the cases of deciduous tree species IF (lai(icir) .GT. 0.1) THEN tree_heights_from_edge(icir) = (28.0/lai(icir)) / mean_height(ipts,ivm,icir) ELSE tree_heights_from_edge(icir) = (28.0/0.1) / mean_height(ipts,ivm,icir) ENDIF ! As suggested by Barry (the father of GALES), the tree height from edge ! should be limited to 9.0 IF (tree_heights_from_edge(icir) .GT. 9.0) tree_heights_from_edge(icir) = 9.0 ! The frame-shaped area around the clear cut gap with X as the width ! of the frame: area_around_gap(icir) = ((length_side_gap + & tree_heights_from_edge(icir)*mean_height(ipts,ivm,icir)*2.0)**2.0) - & clear_cut_max ! Total area of the "frames" in the particular dbh-class of the particula ! PFT:the multiplication by 0.25 is because of the assumption that wind ! direction does not change during the half-hourly time step, therefore ! only one of the four sides of the gaps are affected. ! Use MIN to make sure area_total_closer and area_total_further are ! positive. ! +++JINA This area should not be below zero. I added zero to MIN ! for now, but maybe need to add error to avoid miss ! calculations. area_total_closer(icir) = MIN(n_gap * area_around_gap(icir) * 0.25, & area(ipts)*veget_max(ipts,ivm) - & (n_gap * clear_cut_max), zero) ! The total area which is nor gap neither "frame": area_total_further(icir) = area(ipts) * veget_max(ipts,ivm) - & area_total_closer(icir) - (n_gap * clear_cut_max) IF (area(ipts).GT.zero) THEN ! Get the area fraction closer_area_fraction(icir) = area_total_closer(icir) / (area(ipts)*veget_max(ipts,ivm)) further_area_fraction(icir) = area_total_further(icir) / (area(ipts)*veget_max(ipts,ivm)) ELSE WRITE(numout,*) 'ERROR: the total area of the pixel is zero' CALL ipslerr_p(3,'stomate_windthrow',& 'area of the pixel is zero','','') ENDIF IF(ok_wlsk .AND. ivm .EQ. test_pft .AND. ipts .EQ. test_grid) THEN WRITE(numout,*) 'wlskwlsk_area, ipts, icir',ipts, icir WRITE(numout,*) 'n_gap',n_gap WRITE(numout,*) 'area',area(ipts) WRITE(numout,*) 'area-harvested',area(ipts)-(n_gap*clear_cut_max) WRITE(numout,*) 'gap_area',gap_area_save(ipts,ivm,:) WRITE(numout,*) 'closer_area_f',closer_area_fraction(icir) WRITE(numout,*) 'further_area_f',further_area_fraction(icir) ENDIF !-Debug- IF(printlev_loc .GT. 4 .AND. ivm .EQ. test_pft)THEN WRITE(numout,*) 'Section 2.0: Calculating the gustiness of wind & & in the PFT and area fractions in the grid.' WRITE(numout,*) 'area_closer_fraction, ', closer_area_fraction(icir) WRITE(numout,*) 'area_further_fraction, ', further_area_fraction(icir) ENDIF !- !! 2.2 Calculate ratio to spacing height ! The ratio of spacing to height describes the functional (for ! wind) stand density. If there are many tall trees, the spacing ! is low and thus ratio_spacing_height will be low as well. ! Values of the ratio of mean tree spacing and height are limited ! between 0.075 and 0.45. Bigger values don't make any difference, ! smaller values are not really possible in real life. ! The D/h value are presented in Hale et al., 2015 ratio_spacing_height(icir) = MAX(MIN(current_spacing(icir) / & mean_height(ipts,ivm,icir),0.45),0.075) ! Copy the ratio_spacing_height to the D_H_ratio for output d_h_ratio(ipts,ivm,icir) = ratio_spacing_height(icir) !! 2.3 Calculate functional gap size ! Calculate whether the gap is large enough to increase ! the sensitivity of the remaining trees to wind throw. ! GALES uses the length of the edge of the square-shaped ! gaps calculated in the previous section. gap_size = length_side_gap ! If the gap is wider than 10 times the tree height, ! then wind loading on the stand edge does not increase any ! more. The threshold of 10 is used in the equations below IF (gap_size > 10.0 * mean_height(ipts,ivm,icir)) THEN mean_gap_factor(icir) = 1. max_gap_factor(icir) = 1. ELSE mean_gap_factor(icir) = (0.001 * (gap_size / & mean_height(ipts,ivm,icir))**0.562) / (0.001 * 10.**0.562) max_gap_factor(icir) = (0.0064 * (gap_size / & mean_height(ipts,ivm,icir))**0.3467) / (0.0064 * 10.**0.3467) END IF ! In the following section, there are three types of terms for forest edges: ! the edge (edge), the area around the edge (closer) and the area further ! away from the edge (further). Only closer and further are used in ! the calculation of cws. Gustiness is different in both areas. ! Calculating term_a and term_b for the new gust method. The definition of ! the term_a and term_b are in the Hale et al. (2015) eq(7) ! D/h is the ratio_spacing_height ! x/h is the tree_heights_from_edge term_a(icir) = MAX(-2.1 * ratio_spacing_height(icir) + 0.91, zero) term_b(icir) = 1.0611 * LOG(ratio_spacing_height(icir)) + 4.2 ! New gust factors: ! The gustiness away from the edge is always higher than at the edge. ! We applied term_a*tree_heights_from_edge + term_b for gustiness calculation. ! The gustiness away from the edge is always higher than at the edge. ! For the further area, tree_height_from_edge is always greater than 9.0 ! thus, the above calculation will produce an even higher value for gustiness ! calculation. So, We set tree_heights_from_edge as 9.0. ! The gustiness at the further area will than be higher then the gustiness at ! the closer area. new_gust_closer(icir) = term_a(icir) * tree_heights_from_edge(icir) + term_b(icir) new_gust_further(icir) = term_a(icir) * 9.0 + term_b(icir) new_gust_edge(icir) = term_a(icir) * zero + term_b(icir) ! Calculating gustiness: g_closer(ipts,ivm,icir) = new_gust_closer(icir) g_further(ipts,ivm,icir) = new_gust_further(icir) ! The following section is used to calculate the edge_factor. This section calculates ! how the mean loading on a tree changes as a function of tree height, spacing, distance ! from edge and size of any upwind gap. An upwind gap greater than 10 tree heights ! is assumed to be an infinite gap. The output (edge_factor) is used to modify the ! calculated wind loading on the tree, which assumes a steady wind and a position well ! inside the forest. mean_bm_gust_factor_closer(icir) = (0.68 * ratio_spacing_height(icir) - & 0.0385) + (-0.68 * ratio_spacing_height(icir) + 0.4785) * ((1.7239 * & ratio_spacing_height(icir) + 0.0316)**tree_heights_from_edge(icir)) * & mean_gap_factor(icir) mean_bm_gust_factor_further(icir) = (0.68 * ratio_spacing_height(icir) - & 0.0385) + (-0.68 * ratio_spacing_height(icir) + 0.4785) * ((1.7239 * & ratio_spacing_height(icir) + 0.0316)**9.0) * mean_gap_factor(icir) edge_factor(icir) = mean_bm_gust_factor_closer(icir) / mean_bm_gust_factor_further(icir) ! Copy array for output edge_factor_out(ipts,ivm,icir) = edge_factor(icir) ! Calculating gust_factor accounting for gustiness and edge_factor: gust_factor_closer(ipts,ivm,icir) = g_closer(ipts,ivm,icir) * edge_factor(icir) ! Suggested by Barry ! For further we force edge_factor = 1. This is because everything is ! compared to wind load on tree in the middle of the forest edge_factor(icir) = 1.0 gust_factor_further(ipts,ivm,icir) = g_further(ipts,ivm,icir) * edge_factor(icir) !-Debug- IF(printlev_loc .GT. 4 .AND. ivm .EQ. test_pft)THEN WRITE(numout,*) 'End of Section 2.0: Calculating the & & gustiness of wind in the PFT' WRITE(numout,*) 'term_a, ', term_a(icir) WRITE(numout,*) 'term_b, ', term_b(icir) WRITE(numout,*) 'new_gust_closer, ', new_gust_closer(icir) WRITE(numout,*) 'ratio_spacing_height, ', ratio_spacing_height(:) WRITE(numout,*) 'new_gust_further, ', new_gust_further(icir) WRITE(numout,*) 'new_gust_edge, ', new_gust_edge(icir) WRITE(numout,*) 'tree_heights_from_edge, ', tree_heights_from_edge(icir) WRITE(numout,*) 'Correction factor for the forest & & edge dimension in number of circumference' WRITE(numout,*) 'edge_factor, ', edge_factor(icir) WRITE(numout,*) 'gustness_closer, ', g_closer(ipts,ivm,icir) WRITE(numout,*) 'gustness_futher, ', g_further(ipts,ivm,icir) WRITE(numout,*) 'gust_factor_closer, ' , gust_factor_closer(ipts,ivm,icir) WRITE(numout,*) 'gust_factor_further, ', gust_factor_further(ipts,ivm,icir) WRITE(numout,*) 'Edge factor:', edge_factor(icir) WRITE(numout,*) 'Tree height from edge:', tree_heights_from_edge(icir) WRITE(numout,*) 'Mean bm_gust_factor_closer:', & mean_bm_gust_factor_closer(icir) WRITE(numout,*) 'Mean bm_gust_factor_further:', & mean_bm_gust_factor_further(icir) WRITE(numout,*) 'Mean gap_factor:', mean_gap_factor(icir) ENDIF !- !! 3. Calculating critical moments (tree mechanics section) ! overturning_moment_multiplier: this comes form the species parameters file and ! is related to the depth and type of the soil. Rooting depth has three categories: ! Shallow (less then 0.8 m deep); Deep (deeper than 0.8 m); Average (average values ! to be used when rooting depth is unknown). Soil type has four categories: ! Free-draining mineral soils; Gleyed mineral soils; Peaty mineral soils; Deep peats. ! As these are generic soil types, a modified soil map of Europe would be needed ! for the optimal use of WINDTHROW. IF (soil_type(ipts) == ifree_draining .AND. Rooting_depth(ipts,ivm) == ishallow & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_shallow(ivm) ELSE IF (soil_type(ipts) == ifree_draining .AND. rooting_depth(ipts,ivm) == ishallow & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_shallow_leafless(ivm) ELSE IF (soil_type(ipts) == ifree_draining .AND. rooting_depth(ipts,ivm) == ideep & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR.& plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_deep(ivm) ELSE IF (soil_type(ipts) == ifree_draining .AND. rooting_depth(ipts,ivm) == ideep & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_deep_leafless(ivm) ELSE IF (soil_type(ipts) == ifree_draining .AND. rooting_depth(ipts,ivm) == iaverage & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_average(ivm) ELSE IF (soil_type(ipts) == ifree_draining .AND. rooting_depth(ipts,ivm) == iaverage & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_average_leafless(ivm) ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == ishallow & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_shallow(ivm) ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == ishallow & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_shallow_leafless(ivm) ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == ideep & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_deep(ivm) ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == ideep & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_deep_leafless(ivm) ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == iaverage & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_average(ivm) ELSE IF (soil_type(ipts) == igleyed .AND. rooting_depth(ipts,ivm) == iaverage & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_average_leafless(ivm) ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == ishallow & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peaty_shallow(ivm) ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == ishallow & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peaty_shallow_leafless(ivm) ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == ideep & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peaty_deep(ivm) ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == ideep & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peaty_deep_leafless(ivm) ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == iaverage & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peaty_average(ivm) ELSE IF (soil_type(ipts) == ipeaty .AND. rooting_depth(ipts,ivm) == iaverage & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peaty_average_leafless(ivm) ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == ishallow & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peat_shallow(ivm) ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == ishallow & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peat_shallow_leafless(ivm) ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == ideep & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peat_deep(ivm) ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == ideep & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peat_deep_leafless(ivm) ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == iaverage & .AND. (plant_status(ipts,ivm).EQ.icanopy .OR. & plant_status(ipts,ivm).EQ.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peat_average(ivm) ELSE IF (soil_type(ipts) == ipeat .AND. rooting_depth(ipts,ivm) == iaverage & .AND. (plant_status(ipts,ivm).NE.icanopy .OR. & plant_status(ipts,ivm).NE.ipresenescence)) THEN overturning_moment_multiplier(ipts,ivm) = overturning_peat_average_leafless(ivm) ELSE !average value if parameters are missing overturning_moment_multiplier(ipts,ivm) = 125.0 END IF ! Maximum moment calculated for tree overturning. max_overturning_moment(ipts,ivm,icir) = overturning_moment_multiplier(ipts,ivm) * & stem_mass(ipts,ivm,icir) ! Maximum moment calculated for stem breakage. max_breaking_moment(ipts,ivm,icir) = modulus_rupture(ivm) * f_knot(ivm) * pi * & (mean_dbh(ipts,ivm,icir)**3.0) / 32.0 ! This is where the gustiness of wind is accounted for so that half hourly mean ! wind speeds can be used in the simulations. ! In the end, moments and hence, critical wind speeds are calculated ! for both areas (closer & further) with different gustiness. overturning_moment_closer(ipts,ivm,icir) = max_overturning_moment(ipts,ivm,icir) / & gust_factor_closer(ipts,ivm,icir) overturning_moment_further(ipts,ivm,icir) = max_overturning_moment(ipts,ivm,icir) / & gust_factor_further(ipts,ivm,icir) breaking_moment_closer(ipts,ivm,icir) = max_breaking_moment(ipts,ivm,icir) / & gust_factor_closer(ipts,ivm,icir) breaking_moment_further(ipts,ivm,icir) = max_breaking_moment(ipts,ivm,icir) / & gust_factor_further(ipts,ivm,icir) !-Debug- IF(printlev_loc .GT. 4 .AND. ivm .EQ. test_pft)THEN WRITE(numout,*) 'PFT type, ', ivm WRITE(numout,*) 'POINT id, ', ipts WRITE(numout,*) 'SENSCENCE type,', plant_status(ipts,ivm) WRITE(numout,*) 'Overturning multiplier for free drainage average root leafless all PFTs,', & overturning_free_draining_average_leafless(:) WRITE(numout,*) 'Overturning multiplier for free drainage average root leaf all PFTs,', & overturning_free_draining_average(:) WRITE(numout,*) 'Overturning multiplier for free drainage shallow root leafless all PFTs,', & overturning_free_draining_shallow_leafless(:) WRITE(numout,*) 'Overturning multiplier for free drainage shallow root leaf all PFTs,', & overturning_free_draining_shallow(:) WRITE(numout,*) 'Overturning multiplier for free drainage deep root leafless all PFTs,', & overturning_free_draining_deep_leafless(:) WRITE(numout,*) 'Overturning multiplier for free drainage deep root leaf all PFTs,', & overturning_free_draining_deep(:) WRITE(numout,*) 'End of Section 3 : Calculating critical moments (tree mechanics section)' WRITE(numout,*) 'Diameter class:', icir WRITE(numout,*) 'Overturning moment multiplier, ', & overturning_moment_multiplier(ipts,ivm) WRITE(numout,*) 'maximum moment for overturning, ', max_overturning_moment(ipts,ivm,icir) WRITE(numout,*) 'maximum moment for stem breakag, ', max_breaking_moment(ipts,ivm,icir) WRITE(numout,*) 'overturning closer,', overturning_moment_closer(ipts,ivm,icir) WRITE(numout,*) 'breaking closer, ', breaking_moment_closer(ipts,ivm,icir) WRITE(numout,*) 'overturning further, ', overturning_moment_further(ipts,ivm,icir) WRITE(numout,*) 'breaking further, ', breaking_moment_further(ipts,ivm,icir) ENDIF !- !! 4. Calculating critical wind speed and damage rate !! 4.1 Calculating critical wind speed in the forest and 10 m above the zero-plane displacement ! +++CHECK+++ ! Ideally ORCHIDEE should use only one surface roughness (z0) and displacement ! height. This is not straightforward because here surface roughness and ! displacement height are calculated as a function of wind speed to account for ! streamlining of the canopy under strong winds. This is done in the subroutine ! streamlining(...) in the calculate_wind_process(). Nevertheless, the subroutine ! streamlining uses the canopy structure from ORCHIDEE and when storm damage ! occurs, the canopy structure is adjusted and this new structure will be used ! in condveg.f90 . The inconsistency is thus between the roughness length ! calculations in condveg.f90 and stomate_windthrow.f90. ! +++++++++++ ! for the edge area ! Critical wind speed for stem breakage in the edge area (m/s) distance_from_force = 0.0 call calculate_wind_process(breaking_moment_closer(ipts,ivm,icir), & current_spacing(icir), canopy_depth(icir), & canopy_breadth(icir), canopy_height(icir), & streamlining_n, streamlining_c, & cws_break_closer(ipts,ivm,icir), & z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir),& distance_from_force) ! Critical wind speed for stem breakage 10 m above ! the zero-plane displacement in the edge area(m/s) call elevate(cws_break_closer(ipts,ivm,icir), & z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir), canopy_height(icir), & cws_break_closer_10(ipts,ivm,icir)) ! Critical wind speed for tree overturning in the edge area (m/s) distance_from_force = 0.0 call calculate_wind_process(overturning_moment_closer(ipts,ivm,icir), & current_spacing(icir), canopy_depth(icir), & canopy_breadth(icir), canopy_height(icir), & streamlining_n, streamlining_c, & cws_overturn_closer(ipts,ivm,icir), & z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir),& distance_from_force) ! Critical wind speed for tree overturning 10 m above ! the zero-plane displacement in the edge area (m/s) call elevate(cws_overturn_closer(ipts,ivm,icir), & z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir), canopy_height(icir), & cws_overturn_closer_10(ipts,ivm,icir)) ! IF the soil at 80 below ground is frozen, the module only allows the stem breakage IF ( soil_temp_daily(ipts) .GT. 273.15 ) THEN ! Calculating damage severity in case any of the critical wind speeds is reached ! in the area around gaps (forest edge area). Whichever of the two damage types ! (overturning, stem breakage) has the lowest critical wind speed, it will be the ! occurring one. IF (cws_overturn_closer_10(ipts,ivm,icir) < cws_break_closer_10(ipts,ivm,icir)) THEN cws_final_closer(ipts,ivm,icir) = cws_overturn_closer_10(ipts,ivm,icir) wind_damage_type_closer(ipts,ivm,icir) = ioverturning ELSE cws_final_closer(ipts,ivm,icir) = cws_break_closer_10(ipts,ivm,icir) wind_damage_type_closer(ipts,ivm,icir) = ibreakage END IF ELSE ! when the soil at 80 cm is frozen, the only damage type is limited to the ! stem_breakage cws_final_closer(ipts,ivm,icir) = cws_break_closer_10(ipts,ivm,icir) wind_damage_type_closer(ipts,ivm,icir) = ibreakage ENDIF ! Based on observations of wind damage in Scotland and how much we had to ! adjust the critical wind speed to get agreement between the wind speed ! above the canopy, the calculated critical wind speed and the presence ! of damage (Hale et al, 2015). I am suggesting the relationship below to ! calculate the percentage of damage at a grid point as a function of the ! meteorological wind speed at 10m above the zero-plane displacement ! (prezhd_wind), and the critical wind speed for damage (crit_wind) ! calculated by ForestGALES. I would set maximum to 80% for the moment ! (sets the maximum total amount of damage = 80%) and the scaling factor ! (s) to 0.8 (this affects the rate at which damage percentage increases ! with increasing wind speed). We can adjust these if the seem to be ! wrong or I find other data that I can use to make a better assessment. ! We apply a sigmodial function to tune the "damage level" by using a ! function which is dependent on the difference between the actual wind ! speed and the final_cws (the lower one of the cws from overturning or ! stem_breakage. wind_damage_rate_closer(ipts,ivm,icir) = max_damage_closer(ivm) * ( & (1./(1.+exp(-((max_wind_speed_storm(ipts)- cws_final_closer(ipts,ivm,icir)) / & sfactor_closer(ivm))))) - & (1./(1.+exp(cws_final_closer(ipts,ivm,icir)/sfactor_closer(ivm)))) ) !-Debug- IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN WRITE(numout,*) 'max_damage_closer:', max_damage_closer(ivm) WRITE(numout,*) 'sfactor_closer:', sfactor_closer(ivm) WRITE(numout,*) 'max_wind_speed_storm (m/s):', max_wind_speed_storm(ipts) WRITE(numout,*) 'cws_closer (m/s)', cws_final_closer(ipts,ivm,icir) WRITE(numout,*) 'wind_damage_rate:', wind_damage_rate_closer(ipts,ivm,icir) ENDIF !- !! 4.2 Calculating critical wind speed in the forest and 10 m above the !! zero-plane displacement for the further forest area ! Whichever of the two damage types (overturning, stem breakage) has the ! lowest critical wind speed, it will be the occurring one. ! Critical wind speed for stem breakage in the inner forest area (m/s) distance_from_force = 0.0 call calculate_wind_process(breaking_moment_further(ipts,ivm,icir), & current_spacing(icir), canopy_depth(icir), & canopy_breadth(icir), canopy_height(icir), & streamlining_n, streamlining_c, & cws_break_further(ipts,ivm,icir), & z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir),& distance_from_force) ! Critical wind speed for stem breakage 10 m above the zero-plane ! displacement in the inner forest area (m/s) call elevate(cws_break_further(ipts,ivm,icir), & z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir), canopy_height(icir), & cws_break_further_10(ipts,ivm,icir)) ! Critical wind speed for tree overturning in the inner forest area (m/s) distance_from_force = 0.0 call calculate_wind_process(overturning_moment_further(ipts,ivm,icir),& current_spacing(icir), canopy_depth(icir), & canopy_breadth(icir), canopy_height(icir), & streamlining_n, streamlining_c, & cws_overturn_further(ipts,ivm,icir),& z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir),& distance_from_force) ! Critical wind speed for tree overturning 10 m above the zero-plane ! displacement in the inner forest area (m/s) call elevate(cws_overturn_further(ipts,ivm,icir), & z0_wind(ipts,ivm,icir), zhd_wind(ipts,ivm,icir), canopy_height(icir), & cws_overturn_further_10(ipts,ivm,icir)) ! Whichever of the two damage types (overturning, stem breakage) has the ! lowest critical wind speed, it will be the occurring one. IF (cws_overturn_further_10(ipts,ivm,icir) < cws_break_further_10(ipts,ivm,icir)) THEN cws_final_further(ipts,ivm,icir) = cws_overturn_further_10(ipts,ivm,icir) wind_damage_type_further(ipts,ivm,icir) = ioverturning ELSE cws_final_further(ipts,ivm,icir) = cws_break_further_10(ipts,ivm,icir) wind_damage_type_further(ipts,ivm,icir) = ibreakage END IF ! Wind damage level is considered as a logistic function with sigmoidal ! shape, which is calculated as function of a critical wind speed, an actual ! wind speed, a scaling factor(s)=0.8 and a maximum damage rate(80%)=0.8. wind_damage_rate_further(ipts,ivm,icir) = max_damage_further(ivm) * ( & (1./(1.+exp(-((max_wind_speed_storm(ipts)- cws_final_further(ipts,ivm,icir)) / & sfactor_further(ivm))))) - & (1./(1.+exp(cws_final_further(ipts,ivm,icir)/sfactor_further(ivm)))) ) !-Debug- IF (printlev_loc .GT. 4 .AND. ivm .EQ. test_pft) THEN WRITE(numout,*) 'max_damage_further:', max_damage_further(ivm) WRITE(numout,*) 'sfactor_further:', sfactor_further(ivm) WRITE(numout,*) 'max_wind_speed_storm (m/s):', max_wind_speed_storm(ipts) WRITE(numout,*) 'cws_further (m/s)', cws_final_further(ipts,ivm,icir) WRITE(numout,*) 'wind_damage_rate:', wind_damage_rate_further(ipts,ivm,icir) ENDIF !- !! 5. Determine which trees are killed by wind damage ! We only do the damage for the virtual_circ_class_n GT min_stomate IF (virtual_circ_class_n(ipts,ivm,icir) .GT. min_stomate) THEN !! 5.1. Kill trees in closer area IF (wind_damage_type_closer(ipts,ivm,icir) == ibreakage) THEN kill_break(ipts,ivm,icir) = & wind_damage_rate_closer(ipts,ivm,icir) * circ_class_n(ipts,ivm,icir) * & closer_area_fraction(icir) ! We assume that there is Only one damage type will occur due to ! its lower limiting wind speed kill_uproot(ipts,ivm,icir) = zero ELSE kill_uproot(ipts,ivm,icir) = & wind_damage_rate_closer(ipts,ivm,icir) * & circ_class_n(ipts,ivm,icir) * closer_area_fraction(icir) ! We assume that there is Only one damage type will occur due to ! its lower limited wind speed kill_break(ipts,ivm,icir) = zero ENDIF !! 5.2 Kill trees in further area IF (wind_damage_type_further(ipts,ivm,icir) == ibreakage) THEN ! Total damage = closer + further kill_break(ipts,ivm,icir) = kill_break(ipts,ivm,icir) + & wind_damage_rate_further(ipts,ivm,icir) * & circ_class_n(ipts,ivm,icir) * further_area_fraction(icir) ! We assume that there is only one damage type will occur due to ! its lower limited wind speed kill_uproot(ipts,ivm,icir) = zero ELSE ! Total damage = closer + further kill_uproot(ipts,ivm,icir) = kill_uproot(ipts,ivm,icir) + & wind_damage_rate_further(ipts,ivm,icir) * circ_class_n(ipts,ivm,icir) * & further_area_fraction(icir) ! We assume that there is only one damage type will occur due to ! its lower limited wind speed kill_break(ipts,ivm,icir) = zero ENDIF ELSE ! set no damage for the no vegetation area kill_break(ipts,ivm,icir) = zero kill_uproot(ipts,ivm,icir) = zero END IF ! For the tree height below five meters (7.0m), we assumed that the damage will not ! occur. At the moment, we don't have good solution for the calculation ! of CWS for the small trees, so we just simply keep the small trees. IF ( mean_height(ipts,ivm,icir) .LT. 7.0 ) THEN kill_break(ipts,ivm,icir) = zero kill_uproot(ipts,ivm,icir) = zero ENDIF IF (kill_break(ipts,ivm,icir) .LT. min_stomate) THEN kill_break(ipts,ivm,icir) = zero ENDIF IF (kill_uproot(ipts,ivm,icir) .LT. min_stomate) THEN kill_uproot(ipts,ivm,icir) = zero ENDIF IF(ok_wlsk .AND. ivm .EQ. test_pft .AND. ipts .EQ. test_grid) THEN WRITE(numout,*) 'damage_further',wind_damage_rate_further(ipts,ivm,icir) WRITE(numout,*) 'damage_closer',wind_damage_rate_closer(ipts,ivm,icir) WRITE(numout,*) 'damage_type',wind_damage_type_further(ipts,ivm,icir), & wind_damage_type_closer(ipts,ivm,icir) WRITE(numout,*) 'ccn',circ_class_n(ipts,ivm,icir) ENDIF ! Save mortality to circ_class_cill circ_class_kill(ipts,ivm,icir,forest_managed(ipts,ivm),icut_storm_break) = & kill_break(ipts,ivm,icir) circ_class_kill(ipts,ivm,icir,forest_managed(ipts,ivm),icut_storm_uproot) = & kill_uproot(ipts,ivm,icir) !-Debug- IF (printlev_loc .GT. 4 .AND. ipts==test_grid .AND. ivm==test_pft) THEN WRITE(numout,*) 'End of Section 4: Calculating critical wind speed and damage rate' WRITE(numout,*) 'critical wind speed further, ', cws_final_further(ipts,ivm,icir) WRITE(numout,*) 'critical wind speed closer, ', cws_final_closer(ipts,ivm,icir) WRITE(numout,*) 'wind damage type further, ', wind_damage_type_further(ipts,ivm,icir) WRITE(numout,*) 'wind damage type closer, ', wind_damage_type_closer(ipts,ivm,icir) WRITE(numout,*) 'forest managed type, ', forest_managed(ipts,ivm) WRITE(numout,*) 'circ_class_n(ipts,ivm,icir),', circ_class_n(ipts,ivm,icir) WRITE(numout,*) 'circ_class_kill(ipts,ivm,icir,:,icut_storm_break,', & circ_class_kill(ipts,ivm,icir,forest_managed(ipts,ivm),icut_storm_break) WRITE(numout,*) 'circ_class_kill(ipts,ivm,icir,:,icut_storm_uproot,', & circ_class_kill(ipts,ivm,icir,forest_managed(ipts,ivm),icut_storm_uproot) WRITE(numout,*) 'circ_class_kill, for icut_1:9_:' WRITE(numout,*) circ_class_kill(ipts,ivm,icir,forest_managed(ipts,ivm),1:9) ENDIF !- ! Copy some values for the output files canopy_depth_out(ipts,ivm,icir) = canopy_depth(icir) tree_h_edge_out(ipts,ivm,icir) = tree_heights_from_edge(icir) mean_gap_f_out(ipts,ivm,icir) = mean_gap_factor(icir) gap_size_out(ipts,ivm) = gap_size END DO ! loop ncirc END DO ! loop nvm END DO ! loop npts !! 5. Writing out history variables ! Write the CWS to stomate_history_file DO icir = 1,ncirc !CWS at further area for each diameter class WRITE(var_name,'(A,I3.3)') 'CWS_FURTHER_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & cws_final_further(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc !CWS at closer area for each diameter class WRITE(var_name,'(A,I3.3)') 'CWS_CLOSER_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & cws_final_closer(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc !CWS at closer area in each diameter class for stem breakage WRITE(var_name,'(A,I3.3)') 'CWS_CLOSER_BK_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & cws_break_closer_10(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc !CWS at closer area in each diameter class for overturning WRITE(var_name,'(A,I3.3)') 'CWS_CLOSER_OV_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & cws_overturn_closer_10(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc !CWS at further area in each diameter class for stem break WRITE(var_name,'(A,I3.3)') 'CWS_FURTHER_BK_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & cws_break_further_10(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc !CWS at further area in each diameter class for overturning WRITE(var_name,'(A,I3.3)') 'CWS_FURTHER_OV_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & cws_overturn_further_10(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc !CWS at closer area in each diameter class for stem breakage WRITE(var_name,'(A,I3.3)') 'CWS_TOP_CLO_BK_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & cws_break_closer(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc !CWS at closer area in each diameter class for overturning WRITE(var_name,'(A,I3.3)') 'CWS_TOP_CLO_OV_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & cws_overturn_closer(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc !CWS at further area in each diameter class for stem break WRITE(var_name,'(A,I3.3)') 'CWS_TOP_FUR_BK_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & cws_break_further(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc !CWS at further area in each diameter class for overturning WRITE(var_name,'(A,I3.3)') 'CWS_TOP_FUR_OV_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & cws_overturn_further(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc !D_H_ratio for each diameter class WRITE(var_name,'(A,I3.3)') 'D_H_RATIO_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & d_h_ratio(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! MEAN_HEIGHT for each diameter class WRITE(var_name,'(A,I3.3)') 'MEAN_HEIGHT_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & mean_height(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! MEAN_DBH for each diameter class WRITE(var_name,'(A,I3.3)') 'MEAN_DBH_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & mean_dbh(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! V_IND virtual individual tree numbers for each diameter class WRITE(var_name,'(A,I3.3)') 'V_IND_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & virtual_circ_class_n(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! G_CLOSER gustiness for closer area in each diameter class WRITE(var_name,'(A,I3.3)') 'G_CLOSER_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & gust_factor_closer(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! G_FURTHER gustiness for further area in each diameter class WRITE(var_name,'(A,I3.3)') 'G_FURTHER_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & gust_factor_further(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! G_FURTHER gustiness for further area in each diameter class WRITE(var_name,'(A,I3.3)') 'STEM_MASS_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & stem_mass(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! roughness length WRITE(var_name,'(A,I3.3)') 'WIND_Z0_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & z0_wind(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! displacement height WRITE(var_name,'(A,I3.3)') 'WIND_D_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & zhd_wind(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! canopy depth WRITE(var_name,'(A,I3.3)') 'CANOPY_DEPTH_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & canopy_depth_out(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! tree heights from edge WRITE(var_name,'(A,I3.3)') 'TREE_H_EDGE_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & tree_h_edge_out(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! mean gap factor WRITE(var_name,'(A,I3.3)') 'MEAN_GAP_F_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & mean_gap_f_out(:,:,icir), npts*nvm, horipft_index) ENDDO ! ! gap size WRITE(var_name,'(A)') 'GAP_SIZE' CALL histwrite_p(hist_id_stomate, var_name, itime, & gap_size_out(:,:), npts*nvm, horipft_index) ! DO icir =1,ncirc ! edge factor WRITE(var_name,'(A,I3.3)') 'EDGE_FACTOR_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & edge_factor_out(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! maximum moment for overturning WRITE(var_name,'(A,I3.3)') 'MAX_OV_MOMENT_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & max_overturning_moment(:,:,icir), npts*nvm, horipft_index) ENDDO ! DO icir =1,ncirc ! maximum moment breakage WRITE(var_name,'(A,I3.3)') 'MAX_BK_MOMENT_',icir CALL histwrite_p(hist_id_stomate, var_name, itime, & max_breaking_moment(:,:,icir), npts*nvm, horipft_index) ENDDO ! ! Daily Maximum wind speed CALL histwrite_p(hist_id_stomate, 'DAILY_MAX_WIND', itime, & wind_speed_daily, npts, hori_index) CALL xios_orchidee_send_field("CWS_FURTHER", cws_final_further) CALL xios_orchidee_send_field("CWS_CLOSWER", cws_final_closer) CALL xios_orchidee_send_field("CWS_CLOSER_BK", cws_break_closer_10) CALL xios_orchidee_send_field("CWS_CLOSER_OV", cws_overturn_closer_10) CALL xios_orchidee_send_field("CWS_FURTHER_BK", cws_break_further_10) CALL xios_orchidee_send_field("CWS_FURTHER_OV", cws_overturn_further_10) CALL xios_orchidee_send_field("GAP_AREA",gap_area_save(:,:,1)) CALL xios_orchidee_send_field("DAILY_MAX_WIND", wind_speed_daily) CALL xios_orchidee_send_field("V_IND", virtual_circ_class_n) CALL xios_orchidee_send_field("KILL_STORM_BK", kill_break) CALL xios_orchidee_send_field("KILL_STORM_OV", kill_uproot) END SUBROUTINE wind_damage !! ============================================================================================= !! SUBROUTINE : calculate_wind_process !! !>\BRIEF : This subroutine calculates the breaking and overturning critical wind !! speeds !! !! DESCRIPTION : None !! !! RECENT CHANGE(S): None !! !! RETURN VALUE : !! !! REFERENCES : Gardiner et al. 2010, hale et al. 2012, Raupach 1994. !! !! FLOWCHART : None !! !_ ============================================================================================= SUBROUTINE calculate_wind_process(critical_moment,& current_spacing, canopy_depth, canopy_breadth, canopy_height, & streamlining_n, streamlining_c, & guess_wind_speed, z0_wind, zhd_wind, distance_from_force) IMPLICIT NONE !! 0. Variable and parameter declaration !! 0.1 Input variables REAL, INTENT(in) :: current_spacing REAL, INTENT(in) :: canopy_depth REAL, INTENT(in) :: canopy_breadth REAL, INTENT(in) :: canopy_height REAL, INTENT(in) :: critical_moment REAL, INTENT(in) :: streamlining_n REAL, INTENT(in) :: streamlining_c REAL, INTENT(in) :: distance_from_force !! 0.2 Output variables REAL, INTENT(out) :: guess_wind_speed !! A guess for the critical wind speed REAL, INTENT(out) :: z0_wind REAL, INTENT(out) :: zhd_wind !! 0.3 Modified variables !! 0.4 Local variables INTEGER :: loop_count !! Counting the iteration INTEGER :: max_count=20 !! Max iteration times set to 20 iterations REAL :: delta !! Difference between the guessed !! and the calculated wind speed REAL :: wind_precision !! How close we should !! get with the iterations REAL :: lambda !! REAL :: lambda_capital REAL :: gamma_solved REAL :: psih REAL :: critical_moment_cws REAL, PARAMETER :: air_density = 1.2250 REAL, PARAMETER :: f_crown_mass = 1.136 REAL, PARAMETER :: f_knots = 1.000 !! 0.5 Variables used in calculate_force subroutine and calculate_bending_moment subroutine !! (Both subroutines are not used in the current code) REAL :: gammaSolved REAL :: d REAL :: force_on_tree REAL :: height_of_force REAL :: bending_moment REAL :: porosity REAL :: crit_wind_speed !_ ============================================================================================= !! 1. Subroutine ! This bit belongs to the original form of doing the iterations. For more info, see the notes below. ! guess_wind_speed = 64.0 ! delta = guess_wind_speed / 2.0 guess_wind_speed = 25.0 delta = guess_wind_speed / 2.0 wind_precision = 0.01 loop_count = 0 DO WHILE (delta > wind_precision) CALL streamlining(guess_wind_speed, current_spacing, & canopy_breadth, canopy_depth, canopy_height, & streamlining_n, streamlining_c, gamma_solved, zhd_wind, z0_wind) ! Calculate critical wind speed to compare with wind speeds used for calculating ! streamlining in the iterative calculation of the final critical wind speed) crit_wind_speed = (gamma_solved / current_spacing) * (critical_moment / & (f_crown_mass * air_density * (zhd_wind - distance_from_force) ))**0.5 ! We use the absolute error as the criteria for stopping the iteration ! loop delta = abs( guess_wind_speed - crit_wind_speed ) guess_wind_speed = crit_wind_speed loop_count = loop_count +1 ! Add the maximum iteration times condition for CWS calculation IF (loop_count .GE. max_count) THEN WRITE(numout,*) 'WARANNING from stomate_windthrow.f90: The maximum iteration times for CWS claculation!!' EXIT ENDIF END DO END SUBROUTINE calculate_wind_process !! =========================================================================================== !! SUBROUTINE : streamlining !! !>\BRIEF : This subroutine calculates the streamlining of tree crowns. Streamlining is !! the change of shape of the crowns due to wind. !! !! DESCRIPTION : None !! !! RECENT CHANGE(S): None !! !! RETURN VALUE : Surface roughness and zero-plane displacement among other variables that !! are used in the other subroutines !! REFERENCES : !! Raupach (1994) Simplified expressions for vegetation roughness length and !! zero-plane displacement as functions of canopy height and area index, Boundary-Layer Meteorology !! October 1994, Volume 71, Issue 1, pp 211-216, doi: 10.1007/BF00709229 !! !! FLOWCHART : None !! !_ ============================================================================================ SUBROUTINE streamlining(guess_wind_speed, current_spacing, & canopy_breadth, canopy_depth, a_mean_height, & streamlining_n, streamlining_c, gamma_solved, zhd_wind, z0_wind) IMPLICIT NONE !! 0. Variable and parameter declaration !! 0.1 Input variables REAL, INTENT(in) :: guess_wind_speed !! the guess wind speed (m s{-1}) REAL, INTENT(in) :: current_spacing !! the average spacing between trees (m) REAL, INTENT(in) :: canopy_breadth !! the maximum width of the canopy (m) REAL, INTENT(in) :: canopy_depth !! the length of the live tree crown (m) REAL, INTENT(in) :: a_mean_height !! the mean canopy height (m) REAL, INTENT(in) :: streamlining_n !! streamlining parameter n (-) see Eq. A12 in Hale et al 2015 REAL, INTENT(in) :: streamlining_c !! streamlining parameter c (-) see Eq. A12 in Hale et al 2015 !! 0.2 Output variables REAL, INTENT(out) :: gamma_solved !! a function for solving surface roughness (-) REAL, INTENT(out) :: zhd_wind !! zero plane of displacement height (m) REAL, INTENT(out) :: z0_wind !! surface roughness (m) !! 0.3 Modified variables !! 0.4 Local variables REAL :: lambda !! roughness density, the frontal !! area of roughness elements REAL :: lambda_capital !! canopy area index REAL :: psih !! stability correction function for heat transfer (-) REAL :: porosity !! the effective drag ?? (-) !! see equation A12 in the supplementary material of !! Hale et al (2015) !! The aerodynamic parameter values used here are based on the reference of Raupach (1994) !! --- TEMP ---- !! Need to be discussed for the consistency in the ORCHIDEE !! ------------ REAL, PARAMETER :: c_displacement = 7.5 REAL, PARAMETER :: c_surface = 0.003 REAL, PARAMETER :: c_drag = 0.3 REAL, PARAMETER :: c_roughness = 2.0 REAL, PARAMETER :: lambda_capital_max = 0.6 REAL, PARAMETER :: wind_streamlining_max = 25.0 REAL, PARAMETER :: wind_streamlining_min = 10.0 !_ ============================================================================================= !! 1. Subroutine IF (guess_wind_speed > wind_streamlining_max) THEN ! Above 25 m/s, tree crowns cannot streamline any further. ! The porosity set to its minimum value 2.35*25^-0.51 ~ 0.46 porosity = streamlining_c * wind_streamlining_max**(-streamlining_n) ELSE IF ( guess_wind_speed < wind_streamlining_min ) THEN ! Below 10 m/s, tree crowns do not streamline and porosity set to its ! maximum value 2.35*10^-0.51 ~ 0.73 porosity = streamlining_c * wind_streamlining_min**(-streamlining_n) ELSE ! Wind speed between 10 m/s and 25 m/s, tree crown do the streamlining porosity = streamlining_c * guess_wind_speed**(-streamlining_n) END IF ! lambda = (canopy_breadth/2.0) * canopy_depth * porosity / current_spacing **2.0 ! Here, we directly calculated the mean width of canopy from the value of "canopy_breadth" ! The shape of tree crown is not rhomboid anymore. In the ORCHIDEE-CAN ! the shape of tree crown is circle, so the frontal area should be ! calculated as (pi/4 * canopy_breadth * canopy_depth). lambda = canopy_breadth * canopy_depth * porosity / current_spacing **2.0 lambda_capital = lambda IF (lambda_capital > lambda_capital_max) THEN ! Gamma is needed for the calculation of the roughness length. gamma_solved = 1.0 / ((c_surface + c_drag * lambda_capital_max / 2.0)**0.5) ELSE gamma_solved = 1.0 / ((c_surface + c_drag * lambda_capital / 2.0)**0.5) END IF ! Calculation of the zero-plane displacement: ! see equations from A7 to A11 in the supplementary material of Hale et al. (2015) zhd_wind = a_mean_height * (1.0 - ((1.0 - EXP(-(c_displacement * & lambda_capital)**0.5)) / (c_displacement * lambda_capital)**0.5)) ! Calculation of the aerodynamic roughness length: ! psih = "Profile influence function" at h (height of the roughness elements). psih = LOG(c_roughness) -1.0 + c_roughness**(-1.0) z0_wind = (a_mean_height - zhd_wind) * EXP(psih - ct_karman * gamma_solved) END SUBROUTINE streamlining !! !============================================================================================= !! SUBROUTINE : elevate !! !>\BRIEF : This subroutine converts the critical wind speeds !calculated above to the !! wind speed 10 meters above the zero-plane displacement !! !! DESCRIPTION : None !! !! RECENT CHANGE(S): None !! !! RETURN VALUE : !! !! REFERENCES : !! !! FLOWCHART : None !! !_ !============================================================================================= SUBROUTINE elevate(uh_speed, z0_wind, zhd_wind, mean_height,elevate_result) IMPLICIT NONE !! 0. Variable and parameter declaration !! 0.1 Input variables REAL, INTENT(in) :: uh_speed REAL, INTENT(in) :: z0_wind REAL, INTENT(in) :: zhd_wind REAL, INTENT(in) :: mean_height !! 0.2 Output variables REAL, INTENT(out) :: elevate_result !! 0.3 Modified variables !! 0.4 Local variables !_ !============================================================================================== !! 1. Subroutine elevate_result = uh_speed * log(10.0 / z0_wind) / & log((mean_height - zhd_wind) / z0_wind) END SUBROUTINE elevate !============================================================================================== ! Below subroutines are not used for now !============================================================================================== !! ============================================================================================ !! SUBROUTINE : calculate_force !! !>\BRIEF : This subroutine calculates the force affecting the tree. !! !! DESCRIPTION : None !! !! RECENT CHANGE(S): None !! !! RETURN VALUE : !! !! REFERENCES : !! !! FLOWCHART : None !! !_ ============================================================================================ SUBROUTINE calculate_force(guess_wind_speed, a_gamma_solved, & a_zhd_wind, current_spacing, a_mean_height, & force_on_tree, height_of_force) IMPLICIT NONE !! 0. Variable and parameter declaration !! 0.1 Input variables REAL, INTENT(in) :: guess_wind_speed REAL, INTENT(in) :: a_gamma_solved REAL, INTENT(in) :: a_zhd_wind REAL, INTENT(in) :: current_spacing REAL, INTENT(in) :: a_mean_height !! 0.2 Output variables REAL, INTENT(out) :: force_on_tree REAL, INTENT(out) :: height_of_force !! 0.3 Modified variables !! 0.4 Local variables REAL,PARAMETER :: air_density = 1.2226 !!This should call from constant.f90 (kg m^{-3}) !_ ============================================================================================= !! 1. Subroutine force_on_tree = air_density * (guess_wind_speed * & current_spacing / a_gamma_solved)**2.0 !! --- TEMP --- !! I checked with the reference the force of wind can be calculated as !! air_density*ustar^{2}*D^{2} as the gamma is defined as U/ustar , D is the current spacing !! the force of wind should be as !! force_on_tree = air_density * (U/gamma)^{2} * D^{2} !! ------------ height_of_force = a_zhd_wind ! The maximum of the zero-plane displacement height is locked at 0.8 tree height. ! IF (height_of_force > 0.8 * a_mean_height) THEN ! height_of_force = 0.8 * a_mean_height ! END IF END subroutine calculate_force !! ============================================================================================= !! SUBROUTINE : calculate_bending_moment !! !>\BRIEF : Subroutine to calculate the bending moment !! (to compare with critical moments in the iterative calculation of the critical !! wind speeds) !! !! DESCRIPTION : None !! !! RECENT CHANGE(S): None !! !! RETURN VALUE : !! !! REFERENCES : !! !! FLOWCHART : None !! !_ ============================================================================================= SUBROUTINE calculate_bending_moment (a_force_on_tree, a_height_of_force, bending_moment) IMPLICIT NONE !! 0. Variable and parameter declaration !! 0.1 Input variables REAL, INTENT(in) :: a_force_on_tree REAL, INTENT(in) :: a_height_of_force !! 0.2 Output variables REAL, INTENT(out) :: bending_moment !! 0.3 Modified variables !! 0.4 Local variables REAL :: f_crown_mass !_ ============================================================================================= !! 1. Subroutine bending_moment = a_force_on_tree * a_height_of_force * f_crown_mass END SUBROUTINE calculate_bending_moment END MODULE stomate_windthrow