Context Navigation

stomate_windthrow.f90 @ 7346

Last change on this file since 7346 was 5691, checked in by sebastiaan.luyssaert, 6 years ago

DEV: tested with 13, 37 and 64 PFTs with LCC on different pixels. Some configuration run for 20 years on a given pixel, other crash on another pixel. There is a mass balance problem in sapiens_lcc (ticket #482). This commit fixes a problem with PFT1 in littercalc. This PFT is now fully integrated in LCC and subsequent litter and soil dynamics. veget_max was changed in veget_cov_max where appropriate, a typo in enerbil was corrected.

File size: 109.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: branches/publications/ORCHIDEE_CN_CAN_r5698/src_stomate/stomate_windthrow.f90 @ 7346

Download in other formats: