Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: branches/publications/ORCHIDEE_CAN_r3069/src_stomate/stomate_windfall.f90

Last change on this file was 2966, checked in by sebastiaan.luyssaert, 9 years ago
DEV: add missing windfall module to svn
File size: 76.7 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_wind
3	!
4	! CONTACT : Ferenc.Pasztor@lsce.ipsl.fr
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_windfall.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_windfall
23
24	! modules used:
25	USE constantes
26	USE constantes_soil
27	USE pft_parameters
28	USE function_library
29
30	IMPLICIT NONE
31
32	! private & public routines
33
34	PRIVATE
35	PUBLIC wind_damage
36	!!$ wind_clear, streamlining, calculate_force, calculate_bending_moment, &
37	!!$ calculate_cws, calculate_wind_process, elevate
38
39	LOGICAL, SAVE :: firstcall = .TRUE. !! first call
40
41	CONTAINS
42
43
44	!! ================================================================================================================================
45	!! SUBROUTINE : wind_clear
46	!!
47	!>\BRIEF Set the firstcall flag to .TRUE. and activate initialization
48	!!
49	!_ ================================================================================================================================
50
51	SUBROUTINE wind_clear
52	firstcall = .TRUE.
53	END SUBROUTINE wind_clear
54
55
56	!! ================================================================================================================================
57	!! SUBROUTINE : wind_damage
58	!!
59	!>\BRIEF This subroutine calculates the critical wind speeds for both uprooting and stem breakage
60	!! for each half-hourly period of the simulation. The algorithm follows Barry Gardiner's wind
61	!! damage risk model called GALES (Hale et al. 2015).
62	!!
63	!! DESCRIPTION : The purpose of this module is to describe the actual sensitivity of a given PFT to the wind that is present
64	!! at the same pixel. For each circumference class in each PFT, two critical wind speed values (CWS) are calculated in every
65	!! half-hourly time step by this subroutine, following the approach of differentiating two types of wind damage of trees:
66	!! - uprooting (the whole tree is removed from the soil together with its root plate)
67	!! - stem breakage (the root plate and the bottom part of the stem remains attached to the ground).
68	!! Whichever threshold (i.e. CWS) is reached first for a tree, the according damage type occurs in the model.
69	!!
70	!! The module completes the following tasks:
71	!! 1. Calculates...
72	!!
73	!! RECENT CHANGE(S): None
74	!!
75	!! MAIN OUTPUT VARIABLE(S): :: Critical wind speed for stem breakage 10 m above the zero-plane displacement (m/s): U_Break_10
76	!! Critical wind speed for tree overturn 10 m above the zero-plane displacement (m/s): U_Overturn_10
77	!!
78	!! REFERENCES :
79	!! - Hale, S., Gardiner, B., Nicoll, B., Taylor, P., and Pizzirani, S. (2015). Comparison and validation of three versions
80	!! of a forest wind risk model. Environmental Modelling and Software. In Press.
81	!! ...
82	!!
83	!! FLOWCHART : None
84	!! \streamlining_n
85	!_ ================================================================================================================================
86
87	SUBROUTINE wind_damage (npts, njsc, circ_class_biomass, snow, &
88	veget_max, circ_class_n)
89	!root_dens, senescence, &
90	! circ_class_biomass, circ_class_n, snowdepth, sph, gap_size, snow_depth)
91	! OUTPUT VARIABLES SHOULD BE LISTED HERE, AS WELL.
92
93	IMPLICIT NONE
94
95
96
97	!! 0. Variable declarations
98
99	!! 0.1 Input variables
100
101	INTEGER(i_std), INTENT(in) :: npts !! Domain size @tex $(unitless)$ @endtex
102	INTEGER(i_std), DIMENSION(:), INTENT(in) :: njsc !! Index of the dominant soil textural class in
103	!! the grid cell @tex $(unitless)$ @endtex
104	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass in of an individual tree in each circ class
105	!! @tex $(m^{-2})$ @endtex
106	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of trees within each circumference
107	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
108	!! May sum to less than unity if the pixel has nobio
109	!! area. (unitless, 0-1)
110	REAL(r_std), DIMENSION(:), INTENT(in) :: snow !! Snow mass @tex $ (kg m^{-2})$ @endtex
111
112
113
114	!! 0.2 Output variables
115
116	! OUTPUT VARIABLES TO BE DECLARED HERE
117
118
119	!! 0.3 Modified variables (INTENT(inout))
120
121
122	!! 0.4 Local variables
123
124	INTEGER(i_std) :: ibdl,ipts,ivm,icir,i !! Index for the loops @tex $(unitless)$ @endtex
125	REAL(r_std), DIMENSION(ncirc,60) :: section_height !! The height of the bottom of the 1-m sections
126	!! of the mean tree @tex $(m)$ @endtex. 60 is the
127	!! max number of sections thus max tree height is 60m
128	REAL(r_std), DIMENSION(ncirc,60) :: section_diam !! The diameter of the 1-m sections of the mean
129	!! tree @tex $(m)$ @endtex. max 60 sections.
130	REAL(r_std), DIMENSION(ncirc,60) :: branch_width !! Width of canopy) of the 1-m sections of the mean
131	!! tree @tex $(m)$ @endtex 60 is the max number of
132	!! sections thus max tree height is 60m
133	REAL(r_std), DIMENSION(nbdl,nvm) :: root_dens !! Root density per soil layer (what proportion
134	!! of the root system is found at the given layer)
135	REAL(r_std), DIMENSION(nbdl+1) :: z_soil !! Variable to store depth of the different
136	!! soil layers (m)
137	REAL(r_std), DIMENSION(nvm) :: rpc !! scaling factor for integral
138	REAL(r_std), DIMENSION(npts) :: soil_type !! Soil type (Free-draining mineral soils;
139	!! Gleyed mineral soils; Peaty mineral soils; and
140	!! deep peat soils. See constantes_soil_var.f90
141	!! for the definitions
142	REAL(r_std), DIMENSION(ncirc) :: mean_dbh !! Diameter at breast height (dbh) of the mean tree
143	!! of the circumference-class @tex $(m)$ @endtex
144	REAL(r_std), DIMENSION(ncirc) :: mean_height !! Height of the mean tree of the given dbh-class
145	!! @tex $(m)$ @endtex
146	REAL(r_std), DIMENSION(ncirc) :: lai !! Leaf area index @tex $(unitless)$ @endtex
147	REAL(r_std), DIMENSION(ncirc) :: canopy_breadth !! Canopy breadth (width of canopy) @tex $(m)$ @endtex
148	REAL(r_std), DIMENSION(ncirc) :: canopy_depth !! Canopy depth (height of canopy) @tex $(m)$ @endtex
149	REAL(r_std), DIMENSION(ncirc) :: canopy_height !! The height of the starting point of the canopy from
150	!! below @tex $(m)$ @endtex
151	REAL(r_std), DIMENSION(ncirc) :: mid_canopy !! The middle height of the canopy @tex $(m)$ @endtex
152	REAL(r_std), DIMENSION(ncirc) :: mean_dbh_height !! The height of the diameter (by default it is at
153	!! breast height = 1.3 m) @tex $(m)$ @endtex
154	REAL(r_std), DIMENSION(ncirc) :: snow_mass !! Weight of the snow on the tree @tex $(kg)$ @endtex
155	REAL(r_std), DIMENSION(ncirc) :: diameter_c2 !! Variable for calculating diameters
156	!! @tex $(unitless)$ @endtex
157	REAL(r_std), DIMENSION(ncirc) :: diameter_canopy_base !! Diameter at the base of the canopy
158	!! @tex $(m)$ @endtex
159	REAL(r_std), DIMENSION(ncirc) :: diameter_c1 !! Variable for calculating diameters
160	!! @tex $(unitless)$ @endtex
161	!!$ REAL(r_std), DIMENSION(ncirc) :: no_of_sections !! The number of one-meter sections that the mean tree
162	!!$ !! has @tex $(unitless)$ @endtex
163	REAL(r_std), DIMENSION(ncirc) :: crown_volume !! Crown volume of the mean tree @tex $(m^{3})$ @endtex
164	REAL(r_std), DIMENSION(ncirc) :: crown_mass !! Mass of the mean tree crown @tex $(kg)$ @endtex
165	REAL(r_std), DIMENSION(ncirc) :: stem_mass !! Mass of the mean tree stem @tex $(kg)$ @endtex
166	REAL(r_std), DIMENSION(ncirc) :: current_spacing !! Spacing between trees in the given dbh-class
167	!! @tex $(m)$ @endtex
168	REAL(r_std) :: rooting_depth !! Rooting depth of the mean tree (three categories;
169	!! 1: shallow, 2: deep, 3:average
170	!! See constantes_soil_var.f90 for the definitions
171	!!$ REAL(r_std) :: canopy_density !! Density of the tree canopy, i.e. of the
172	!!$ !! branches and leaves @tex $(kg/m^{3})$ @endtex
173	REAL(r_std) :: streamlining_c !! Streamlining parameter. @tex $(unitless)$ @endtex.
174	!! Streamlining is the change of shape of the crowns
175	!! due to wind.
176	REAL(r_std) :: streamlining_n !! Streamlining parameter. @tex $(unitless)$ @endtex.
177	!! Streamlining is the change of shape of the crowns
178	!! due to wind.
179	REAL(r_std) :: streamlining_rb !! Streamlining parameter. @tex $(unitless)$ @endtex.
180	!! Streamlining is the change of shape of the crowns
181	!! due to wind.
182
183
184
185
186
187	REAL(r_std), DIMENSION(ncirc) :: area_timber_removals_5_years !! Area of the total timber removal in the
188	!! previous 5 year @tex $(m^{2})$ @endtex
189	REAL(r_std), DIMENSION(ncirc) :: area_icir !! The total area of the circumference class @tex $(m^{2})$ @endtex
190	REAL(r_std), DIMENSION(ncirc) :: tree_heights_from_edge !! Width of the forest edge area, i.e. area_closer (in number of tree heights)
191	!! @tex $(unitless)$ @endtex
192	REAL(r_std), DIMENSION(ncirc) :: n_gap !! The number of gaps with an area of clear_cut_max @tex $(unitless)$ @endtex
193	REAL(r_std), DIMENSION(ncirc) :: length_side_gap !! The length of one side of the square-shaped gap @tex $(m)$ @endtex
194	REAL(r_std), DIMENSION(ncirc) :: area_around_gap !! The framing edge area around the gap with different gustiness
195	!! @tex $(m^{2})$ @endtex
196	REAL(r_std), DIMENSION(ncirc) :: area_total_closer !! The total area of forest edge areas affected by wind @tex $(m^{2})$ @endtex
197	REAL(r_std), DIMENSION(ncirc) :: area_total_further !! The total area of forest without edge areas @tex $(m^{2})$ @endtex
198	REAL(r_std), DIMENSION(ncirc) :: ratio_spacing_height !! Ratio of tree height and tree spacing @tex $(unitless)$ @endtex
199	REAL(r_std), DIMENSION(ncirc) :: gap_size !! The length (width of the cross-section) of the gap @tex $(m)$ @endtex
200	REAL(r_std), DIMENSION(ncirc) :: ratio_gap_height !! Ratio of gap size and tree height @tex $(unitless)$ @endtex
201	REAL(r_std), DIMENSION(ncirc) :: mean_gap10 !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
202	REAL(r_std), DIMENSION(ncirc) :: max_gap10 !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
203	REAL(r_std), DIMENSION(ncirc) :: mean_gap_factor !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
204	REAL(r_std), DIMENSION(ncirc) :: max_gap_factor !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
205	REAL(r_std), DIMENSION(ncirc) :: mean_bm_closer !! Mean bending moment for the forest edge area @tex $(Nm/kg)$ @endtex
206	REAL(r_std), DIMENSION(ncirc) :: max_bm_closer !! Maximum bending moment for the forest edge area @tex $(Nm/kg)$ @endtex
207	REAL(r_std), DIMENSION(ncirc) :: mean_bm_further !! Mean bending moment for the forest area @tex $(Nm/kg)$ @endtex
208	REAL(r_std), DIMENSION(ncirc) :: max_bm_further !! Maximum bending moment for the forest area @tex $(Nm/kg)$ @endtex
209	REAL(r_std), DIMENSION(ncirc) :: mean_bm_edge !! Mean bending moment for the forest edge @tex $(Nm/kg)$ @endtex
210	REAL(r_std), DIMENSION(ncirc) :: max_bm_edge !! Maximum bending moment for the forest edge @tex $(Nm/kg)$ @endtex
211	REAL(r_std), DIMENSION(ncirc) :: mean_bm_edge_gap !! Mean bending moment for the forest edge - gap factor considered
212	!! @tex $(Nm/kg)$ @endtex
213	REAL(r_std), DIMENSION(ncirc) :: max_bm_edge_gap !! Maximum bending moment for the forest edge - gap factor considered
214	!! @tex $(Nm/kg)$ @endtex
215	REAL(r_std), DIMENSION(ncirc) :: term_a !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
216	REAL(r_std), DIMENSION(ncirc) :: term_b !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
217	REAL(r_std), DIMENSION(ncirc) :: old_gust_closer !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
218	REAL(r_std), DIMENSION(ncirc) :: old_gust_further !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
219	REAL(r_std), DIMENSION(ncirc) :: old_gust_edge !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
220	REAL(r_std), DIMENSION(ncirc) :: old_gust_edge_gap !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
221	REAL(r_std), DIMENSION(ncirc) :: new_gust_closer !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
222	REAL(r_std), DIMENSION(ncirc) :: new_gust_further !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
223	REAL(r_std), DIMENSION(ncirc) :: new_gust_edge !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
224	REAL(r_std), DIMENSION(ncirc) :: new_gust_edge_gap_closer !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
225	REAL(r_std), DIMENSION(ncirc) :: new_gust_edge_gap_further !! Variable used for calculating gustiness @tex $(unitless)$ @endtex
226	REAL(r_std), DIMENSION(ncirc) :: mean_bm_gust_factor_closer !! Mean bending moment of the forest edge area (for calculating
227	!! the edge factor) @tex $(Nm/kg)$ @endtex
228	REAL(r_std), DIMENSION(ncirc) :: mean_bm_gust_factor_further !! Mean bending moment of the inner forest area (for calculating
229	!! the edge factor) @tex $(Nm/kg)$ @endtex
230	REAL(r_std), DIMENSION(ncirc) :: edge_factor !! The ratio of the mean bending moment in the forest edge area
231	!! and mean bending moment in the inner forest area
232	!! @tex $(unitless)$ @endtex
233	REAL(r_std), DIMENSION(ncirc) :: gust_factor_closer !! Gust factor for calculating the critical wind speed of the forest
234	!! edge area @tex $(unitless)$ @endtex
235	REAL(r_std), DIMENSION(ncirc) :: gust_factor_further !! Gust factor for calculating the critical wind speed of the inner
236	!! forest area @tex $(unitless)$ @endtex
237	REAL(r_std), DIMENSION(ncirc) :: g_closer !! Gustiness of the forest edge area @tex $(unitless)$ @endtex
238	REAL(r_std), DIMENSION(ncirc) :: g_further !! Gustiness of the inner forest area @tex $(unitless)$ @endtex
239	REAL(r_std), DIMENSION(nvm) :: overturning_moment_multiplier!! Overturning moment multiplier representing the effects of species,
240	!! soil type, soil depths and senescence @tex $(Nm/kg)$ @endtex
241	REAL(r_std), DIMENSION(ncirc) :: max_overturning_moment !! The maximum overturning moment the mean tree can withstand
242	!! @tex $(Nm/kg)$ @endtex
243	REAL(r_std), DIMENSION(ncirc) :: max_breaking_moment !! The maximum breaking moment the mean tree can withstand
244	!! @tex $(Nm/kg)$ @endtex
245	REAL(r_std), DIMENSION(ncirc) :: overturning_moment_closer !! The maximum overturning moment the mean tree can withstand in the
246	!! forest edge area with gustiness taken into account @tex $(Nm/kg)$ @endtex
247	REAL(r_std), DIMENSION(ncirc) :: overturning_moment_further !! The maximum overturning moment the mean tree can withstand in the
248	!! inner forest area with gustiness taken into account @tex $(Nm/kg)$ @endtex
249	REAL(r_std), DIMENSION(ncirc) :: breaking_moment_closer !! The maximum breaking moment the mean tree can withstand in the
250	!! forest edge area with gustiness taken into account @tex $(Nm/kg)$ @endtex
251	REAL(r_std), DIMENSION(ncirc) :: breaking_moment_further !! The maximum breaking moment the mean tree can withstand in the
252	!! inner forest area with gustiness taken into account @tex $(Nm/kg)$ @endtex
253	REAL(r_std), DIMENSION(ncirc) :: cws_break_closer !! Critical wind speed for stem breakage in the forest edge area
254	!! @tex $(m/s; mean hourly wind speed)$ @endtex
255	REAL(r_std), DIMENSION(ncirc) :: cws_break_closer_10 !! Critical wind speed for stem breakage in the forest edge area
256	!! @tex $(m/s; mean hourly wind speed)$ @endtex
257	REAL(r_std), DIMENSION(ncirc) :: cws_overturn_closer !! Critical wind speed for overturning in the forest edge area
258	!! @tex $(m/s; mean hourly wind speed)$ @endtex
259	REAL(r_std), DIMENSION(ncirc) :: cws_overturn_closer_10 !! Critical wind speed for overturning in the forest edge area
260	!! @tex $(m/s; mean hourly wind speed)$ @endtex
261	REAL(r_std), DIMENSION(ncirc) :: cws_final_closer !! The critical wind speed of the final form of damage (breakage / uprooting)
262	!! in the forest edge area @tex $(m/s; mean hourly wind speed)$ @endtex
263	REAL(r_std), DIMENSION(npts) :: wind_speed_actual !! The actual wind speed over the pixel
264	!! @tex $(m/s; mean hourly wind speed)$ @endtex
265	LOGICAL :: wind_damage_type_uprooting_closer !! The final form of damage (breakage / uprooting) in the forest edge area
266	REAL(r_std), DIMENSION(ncirc) :: wind_damage_rate_closer !! Damage rate (%) in the forest edge area in case of damage
267	!! @tex $(0-1)$ @endtex
268	REAL(r_std), DIMENSION(ncirc) :: cws_break_further !! Critical wind speed for stem breakage in the inner forest area
269	!! @tex $(m/s; mean hourly wind speed)$ @endtex
270	REAL(r_std), DIMENSION(ncirc) :: cws_break_further_10 !! Critical wind speed for stem breakage in the inner forest area
271	!! @tex $(m/s; mean hourly wind speed)$ @endtex
272	REAL(r_std), DIMENSION(ncirc) :: cws_overturn_further !! Critical wind speed for overturning in the inner forest area
273	!! @tex $(m/s; mean hourly wind speed)$ @endtex
274	REAL(r_std), DIMENSION(ncirc) :: cws_overturn_further_10 !! Critical wind speed for overturning in the inner forest area
275	!! @tex $(m/s; mean hourly wind speed)$ @endtex
276	REAL(r_std), DIMENSION(ncirc) :: cws_final_further !! The critical wind speed of the final form of damage (breakage / uprooting)
277	!! in the inner forest area @tex $(m/s; mean hourly wind speed)$ @endtex
278	LOGICAL :: wind_damage_type_uprooting_further !! The final form of damage (breakage / uprooting) in the inner forest area
279	REAL(r_std), DIMENSION(ncirc) :: wind_damage_rate_further !! Damage rate (%) in the forest edge area in case
280	!! of damage in the inner forest area @tex $(0-1)$ @endtex
281
282	!_ ============================================================================================================================================================================
283
284	!+++CHECK+++
285	! We could split this routine and calculate the sections
286	! only once a day and then calculate the critical wind
287	! speeds every half hour. If if the co
288	!+++++++++++
289
290	!! 1. Getting the input variables in the required form for the actual sub-routine
291
292	! Looping over pixels
293	DO ipts = 1, npts
294
295	!! 1.1 Calculate root profile
296	z_soil(1) = 0.
297	z_soil(2:nbdl+1) = diaglev(1:nbdl)
298
299	! Scaling factor for integration
300	rpc(:) = un / ( un - EXP( -z_soil(nbdl) / humcste(:) ) )
301
302	! Calculate root density per layer per pft
303	! Loop over # soil layers
304	DO ibdl = 1, nbdl
305
306	root_dens(ibdl,:) = rpc(:) * &
307	( EXP( -z_soil(ibdl)/humcste(:)) - &
308	EXP( -z_soil(ibdl+1)/humcste(:)) )
309
310	ENDDO ! Loop over # soil layers
311
312	!! 1.2 Defining the soil type
313	! This sub-routine uses 4 types of soil for calculating the
314	! critical wind speed, so 12 different soil types present in
315	! other modules of ORCHIDEE need to be grouped into these
316	! 4 generic types.
317
318	! +++CHECK+++
319	! For the moment only 3 soil types are distinguished.
320	! Needs to be checked with all 12 USDA soil types.
321	! NOTE: your initial approach to use text is correct
322	! but makes it more difficult to use IF statements.
323	! I implemented the classic work around. Defined in
324	! constantes_soil_var.f90. Update the documentation
325	! in that routine.
326	SELECT CASE (njsc(ipts))
327	CASE (1:3)
328	soil_type(ipts) = ifree_draining
329	CASE (4:6)
330	soil_type(ipts) = igleyed
331	CASE (7:9)
332	soil_type(ipts) = ipeaty
333	CASE (10:12)
334	soil_type(ipts) = ipeat
335	CASE DEFAULT
336	CALL ipslerr_p(3, 'stomate_windfall.f90',&
337	'soil type is not defined','','')
338	END SELECT
339	!++++++++++++
340
341	! Looping over PFT-s
342	DO ivm = 2, nvm
343
344	! Don't do the calculations if there is no vegetation
345	! or the vegetation is bare soil, grass or crop
346	IF ((veget_max(ipts,ivm).LT.min_stomate) .OR. .NOT.is_tree(ivm)) THEN
347	CYCLE
348	ENDIF
349
350	!! 1.3 Calculating the rooting depth
351	! If root density (proportion of root mass) is higher
352	! than 3% in the 11th soil layer (1.024-2.056 m) below
353	! the soil surface, then we treat rooting depth as
354	! "deep". Otherwise, it's "shallow". Note that the
355	! threshold of 3% is arbitrary and it might be needed to
356	! justify at some point in the future.
357
358	! +++CHECK+++
359	! It was hard coded (index = 11). We don't like hard coded
360	! even if it is unlikely that the 11-layers will change
361	! any time soon. The depth of the 11th layer depends on
362	! the total depth of the soil (dpu_max). Would be better
363	! to use ::dpu_max rather than assuming the 2 m.
364	! Also, there is new scheme being tested where nbdl may
365	! have changed. CHECK when using the new scheme.
366	IF (root_dens(nbdl,ivm) > 0.03) THEN
367
368	rooting_depth = ideep
369
370	ELSE
371
372	rooting_depth = ishallow
373
374	END IF
375
376	!! 1.4 Set crown parameters
377	! To set crown parameters leaf mass needs to be available
378	IF (SUM(circ_class_biomass(ipts,ivm,:,ileaf,icarbon)) .GT. min_stomate) THEN
379
380	! There is a canopy
381	! +++CHECK+++
382	! We can calculate the canopy density internally. We
383	! have leaf mass and crown volume. Adding part of the
384	! branches would be more challenging but possible
385	!!$ canopy_density = canopy_density_leaf(ivm)
386	streamlining_c = streamlining_c_leaf(ivm)
387	streamlining_n = streamlining_n_leaf(ivm)
388	streamlining_rb = streamlining_rb_leaf(ivm)
389	! +++++++++++
390
391	ELSE
392
393	! If the tree is leafless, the following parameters
394	! become different. The weight of leaves will be absent.
395	! +++CHECK+++
396	! the definition says this includes leaves and branches.
397	! If set to zero (as is the case here) there are no
398	! branches!
399	!!$ canopy_density = canopy_density_leafless(ivm)
400	! +++++++++++
401
402	! Streamlining changes as the surface of leaves is absent.
403	streamlining_c = streamlining_c_leafless(ivm)
404	streamlining_n = streamlining_n_leafless(ivm)
405	streamlining_rb = streamlining_rb_leafless(ivm)
406
407	END IF
408
409	!---TEMP---
410	IF(ld_windfall)THEN
411	WRITE(numout,*) 'Initialisation of critical windspeed'
412	WRITE(numout,*) 'ipts, ivm, ',ipts, ivm
413	WRITE(numout,*) 'root_dens(ibdl,ivm), ',root_dens(:,ivm)
414	WRITE(numout,*) 'soil_type(ipts), ',soil_type(ipts)
415	WRITE(numout,*) 'rooting_depth, ',rooting_depth
416	WRITE(numout,*) 'streamlining_c, ', streamlining_c
417	WRITE(numout,*) 'streamlining_n, ', streamlining_n
418	WRITE(numout,*) 'streamlining_rb, ', streamlining_rb
419	ENDIF
420	!----------
421
422	!! 1.5 Calculate stand characteristics
423	! Calculating diameter at breast height (dbh) in m of the
424	! mean tree of all circ classes.
425	mean_dbh(:) = wood_to_dia(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
426
427	! Calculating the height in m of the mean tree of the given
428	! dbh-class using
429	mean_height(:) = wood_to_height(&
430	circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
431
432	! Stem mass is calculated as the aboveground woody biomass
433	! minus the branches. circ_class_biomass is the biomass of
434	! an individual tree (gC tree-1)
435	stem_mass(:) = (circ_class_biomass(ipts,ivm,:,isapabove,icarbon) + &
436	circ_class_biomass(ipts,ivm,:,iheartabove,icarbon)) * &
437	(un - branch_ratio(ivm))
438
439	! Crown_mass is calculated as the leaf mass and the branch mass)
440	crown_mass(:) = circ_class_biomass(ipts,ivm,:,ileaf,icarbon) + &
441	(circ_class_biomass(ipts,ivm,:,isapabove,icarbon) + &
442	circ_class_biomass(ipts,ivm,:,iheartabove,icarbon)) * &
443	branch_ratio(ivm)
444
445
446	DO icir = ncirc,1,-1
447
448	!+++CHECK+++
449	! add a while loop. if CWS is not exceed for the
450	! largest diameter class then there is no need
451	! to repeat these calculations for the other
452	! dbh classes
453	!+++++++++++
454
455	! Calculating leaf area index
456	lai(icir) = biomass_to_lai(&
457	circ_class_biomass(ipts,ivm,icir,ileaf,icarbon),ivm)
458
459	! Calculating the current tree spacing in the PFT
460	! The current tree spacing between tree stems are calculated
461	! from the number of stems per hectare. circ_class_n is the
462	! number of stems per m2; it's converted to stems per hectare
463	! with m2_to_ha which has a value of 10000.
464	! +++CHECK+++
465	! this does not account for competiton/overlap
466	! expecially for the largest dbh classes the
467	! spacing will be very large. Calculating a
468	! single current spacing at the stand level
469	! probably makes more sense
470	IF (circ_class_n(ipts,ivm,icir) .GT. min_stomate) THEN
471	current_spacing(icir) = 100.0 / &
472	SQRT(m2_to_ha * circ_class_n(ipts,ivm,icir))
473	ELSE
474	current_spacing(icir) = zero
475	ENDIF
476	! +++++++++++
477
478	ENDDO
479
480	! At this point, we have the weight of the stem and the crown, which
481	! is used in calculating the resistance force of the tree against wind.
482	!---TEMP---
483	IF(ld_windfall)THEN
484	WRITE(numout,*) 'Calculated stem sections'
485	WRITE(numout,*) 'mean dbh, ', mean_dbh(:)
486	WRITE(numout,*) 'mean height, ', mean_height(:)
487	WRITE(numout,*) 'stem_mass, ', stem_mass(:)
488	WRITE(numout,*) 'crown_volume, ',crown_volume(:)
489	WRITE(numout,*) 'crown_mass, ',crown_mass(:)
490	WRITE(numout,*) 'lai, ',lai(:)
491	WRITE(numout,*) 'current_spacing, ',current_spacing(:)
492	ENDIF
493	!----------
494
495
496
497
498	!! 3. Calculating the gustiness of wind in the PFT
499
500	!!$ ! IMPORTANT: AS WE DO NOT HAVE INFORMATION ON STAND EDGES IN THE
501	!!$ ! ORCHIDEE SIMULATIONS, A POSSIBLE PROXY TO USE HERE IS A VARIABLE
502	!!$ ! THAT REPRESENTS EARLIER HARVESTS IN THE PREVIOUS 5 YEARS. AFTER
503	!!$ ! SUCH A PERIOD, THE FOREST EDGE IS ASSUMED TO GET ACCLIMATIZED
504	!!$ ! TO THE NEW CIRCUMSTANCES, AND IT'S VULNERABILITY DIMINISHES
505	!!$ ! (see Persson, 1975; Valinger and Fridman, 2011)
506	!!$
507	!!$ ! There is a big increase in gustiness after a certain distance from
508	!!$ ! the edge of the forest, hence we divide the area of the PFT into two
509	!!$ ! smaller areas:
510	!!$ ! - the area which is further than X tree heights from any edge that
511	!!$ ! was created in the previous five years (area_total_further)
512	!!$ ! - the area which is closer than X tree heights from any edge that
513	!!$ ! was created in the previous five years (area_total_closer).
514	!!$
515	!!$ ! Of course, these smaller areas depend on the fragmentation of the
516	!!$ ! gaps (created in the last five years). As we have no information
517	!!$ ! on this, we assume a gap size for each country depending on the
518	!!$ ! maximum allowed size of clear cuts (clear_cut_max) and that the
519	!!$ ! gaps are square-shaped. Then, we divide the area of harvest in the
520	!!$ ! previous five years (area_timber_removals_5_years) by clear_cut_max.
521	!!$ ! Hence, from the number of gaps we can calculate the area_total_closer.
522	!!$ ! As each circumference class will have its forest area divided, we
523	!!$ ! calculate two gustiness values and hence, two critical wind speeds.
524	!!$
525	!!$ ! NOTE: We have to add a forest management map that contains the maximum
526	!!$ ! allowed areas of clear fellings. As this is not available yet, we use
527	!!$ ! an area of 2 ha as gap size for all calculations.
528	!!$
529	!!$ ! The value of X depends on the leaf area index of the PFT and is calculated as follows:
530	!!$ tree_heights_from_edge(:) = 28.0 * mean_height(icir) /&
531	!!$ lai(icir)
532	!!$ ! This comes from the equation of calculating canopy penetration depth in Pan et al. 2015.
533	!!$
534	!!$ ! Calculating the total area of timber removals in the previous five years:
535	!!$ ! NOTE: Could you help with this?
536	!!$ ! +++CHECK+++
537	!!$ ! needs to become a SAVE in stomate. Just gave it a value so that we can compile
538	!!$ area_timber_removals_5_years(icir) = 1.0
539	!!$ ! +++++++++++
540	!!$
541	!!$ ! Calculating the total forest area belonging to the circumference class:
542	!!$ ! NOTE: Could you help with this?
543	!!$ !+++CHECK+++
544	!!$ ! We already discussed this by mail. Area per circ_class is not easy
545	!!$ ! because the crowns are allowed to overlap. We can find a solution.
546	!!$ ! Area of the whole age class is given by
547	!!$ ! area = resolution_x(ipts) * resolution_y(ipts) * veget_max(ipts,ivm)
548	!!$ ! Need to discuss why this is needed. Just want to get the code
549	!!$ ! compiling for now so I gave it a neutral value
550	!!$ area_icir(icir) = 1.0
551	!!$ !++++++++++++
552	!!$
553	!!$ ! Number of gaps with the area of the maximum allowed clearfelling:
554	!!$ n_gap(icir) = INT(area_timber_removals_5_years(icir) / clear_cut_max)
555	!!$ ! The length of the sides of the squared-shaped gaps:
556	!!$ length_side_gap(icir) = SQRT(clear_cut_max)
557	!!$ ! The frame-shaped area around a 2-ha gap with X as the width of the frame:
558	!!$ area_around_gap(icir) = ((length_side_gap(icir) +&
559	!!$ tree_heights_from_edge(icir) * mean_height(icir))**2.0) - clear_cut_max
560	!!$ ! Total area of the "frames" in the particular dbh-class of the particular PFT:
561	!!$ ! The multiplication by 0.25 is because of the assumption that wind direction does not
562	!!$ ! change during the half-hourly time step, therefore only one of the four sides of the
563	!!$ ! gaps are affected.
564	!!$ area_total_closer(icir) = n_gap(icir) * area_around_gap(icir) * 0.25
565	!!$ ! The total area which is nor gap neither "frame":
566	!!$ area_total_further(icir) = area_icir(icir) - area_total_closer(icir) -&
567	!!$ (n_gap(icir) * clear_cut_max)
568	!!$
569	!!$
570	!!$ ! At this point, we know how areas with different gustiness will look like,
571	!!$ ! now we calculate gustiness for these areas.
572	!!$
573	!!$ ! Values of the ratio of mean tree spacing and height are limited
574	!!$ ! between 0.075 and 0.45. Bigger values don't make any difference,
575	!!$ ! smaller values are not really possible in real life.
576	!!$ IF (current_spacing(icir) / mean_height(icir) < 0.075) THEN
577	!!$ ratio_spacing_height(icir) = 0.075
578	!!$ ELSE
579	!!$ IF (current_spacing(icir) / mean_height(icir) > 0.45) THEN
580	!!$ ratio_spacing_height(icir) = 0.45
581	!!$ ELSE
582	!!$ ratio_spacing_height(icir) = current_spacing(icir) / &
583	!!$ mean_height(icir)
584	!!$ END IF
585	!!$ END IF
586	!!$
587	!!$
588	!!$ ! Calculating gap characteristics:
589	!!$
590	!!$ ! An originally input variable in GALES. We use the length of the edge of the square-shaped
591	!!$ ! gaps calculated in the previous section.
592	!!$ gap_size(icir) = length_side_gap(icir)
593	!!$
594	!!$
595	!!$ ! If the gap is wider than 10 times the tree height,
596	!!$ ! then wind loading on the stand edge does not increase any more.
597	!!$ IF (gap_size(icir) > 10.0 * mean_height(icir)) THEN
598	!!$ gap_size(icir) = 10.0 * mean_height(icir)
599	!!$ ELSE
600	!!$ gap_size(icir) = gap_size(icir)
601	!!$ END IF
602	!!$
603	!!$ ! Calculating factors for calculating bending moments:
604	!!$ ratio_gap_height(icir) = gap_size(icir) / mean_height(icir)
605	!!$ mean_gap10(icir) = 0.001 * 10.0**0.562
606	!!$ max_gap10(icir) = 0.0064 * 10.0**0.3467
607	!!$ mean_gap_factor(icir) = (0.001 * ratio_gap_height(icir)**0.562) / &
608	!!$ mean_gap10(icir)
609	!!$ max_gap_factor(icir) = (0.0064 * ratio_gap_height(icir)**0.3467) / &
610	!!$ max_gap10(icir)
611	!!$
612	!!$ ! In the following section, there are two types of terms for forest edges:
613	!!$ ! the edge and the edge area.
614	!!$ ! The first one represents the very edge of the forest (0 m from the gap),
615	!!$ ! the second one represents the area around the gap (area_around_gap) where
616	!!$ ! gustiness is different than in the rest of the forest.
617	!!$
618	!!$ ! In the original code of GALES, there are two methods to calculate gustiness inside
619	!!$ ! the forest stand and they are both used, depending on the values of a few particular
620	!!$ ! variables. That's the reason for the existence of variables e.g. old_gust_edge
621	!!$ ! and new_gust_edge.
622	!!$
623	!!$ ! The next bit corrects the gustiness for the effect of an upwind gap. First, it calculates
624	!!$ ! the gustiness well inside the forest and at the edge, with a 10-tree gap size and with
625	!!$ ! the specified gap size. It does that using the old method of calculating the gustiness
626	!!$ ! because the gap factors are designed for this method. It then uses these to correct the
627	!!$ ! gustiness at the edge calculated using the new method to account for a gap. Gustiness
628	!!$ ! inside the forest is proportionally adjusted to reflect how far back you are (either X or
629	!!$ ! 9 tree heights).
630	!!$ ! At the edge of the forest, with a very large gap you should have less gustiness than well
631	!!$ ! inside the forest.
632	!!$ ! Gustiness is calculated for both the closer and further areas from the gaps,
633	!!$ ! the critical wind speed will differ accordingly.
634	!!$ ! Zeroes in the power term are only used for completeness.
635	!!$ mean_bm_closer(icir) = (0.68 * ratio_spacing_height(icir) - 0.0385) + &
636	!!$ (-0.68 * ratio_spacing_height(icir) + 0.4785) * (1.7239 * &
637	!!$ ratio_spacing_height(icir) + 0.0316)**tree_heights_from_edge(icir)
638	!!$
639	!!$ max_bm_closer(icir) = (2.7193 * ratio_spacing_height(icir) - 0.061) + &
640	!!$ (-1.273 * ratio_spacing_height(icir) + 0.9701) * (1.1127 * &
641	!!$ ratio_spacing_height(icir) + 0.0311)**tree_heights_from_edge(icir)
642	!!$
643	!!$ mean_bm_further(icir) = (0.68 * ratio_spacing_height(icir) - 0.0385) + &
644	!!$ (-0.68 * ratio_spacing_height(icir) + 0.4785) * (1.7239 * &
645	!!$ ratio_spacing_height(icir) + 0.0316)**9.0
646	!!$
647	!!$ max_bm_further(icir) = (2.7193 * ratio_spacing_height(icir) - 0.061) + &
648	!!$ (-1.273 * ratio_spacing_height(icir) + 0.9701) * (1.1127 * &
649	!!$ ratio_spacing_height(icir) + 0.0311)**9.0
650	!!$
651	!!$ mean_bm_edge(icir) = (0.68 * ratio_spacing_height(icir) - 0.0385) + &
652	!!$ (-0.68 * ratio_spacing_height(icir) + 0.4785) * (1.7239 * &
653	!!$ ratio_spacing_height(icir) + 0.0316)**0.0
654	!!$
655	!!$ max_bm_edge(icir) = (2.7193 * ratio_spacing_height(icir) - 0.061) + &
656	!!$ (-1.273 * ratio_spacing_height(icir) + 0.9701) * (1.1127 * &
657	!!$ ratio_spacing_height(icir) + 0.0311)**0.0
658	!!$
659	!!$ mean_bm_edge_gap(icir) = (0.68 * ratio_spacing_height(icir) - 0.0385) + &
660	!!$ (-0.68 * ratio_spacing_height(icir) + 0.4785) * ((1.7239 * &
661	!!$ ratio_spacing_height(icir) + 0.0316)*0.0) mean_gap_factor(icir)
662	!!$
663	!!$ max_bm_edge_gap(icir) = (2.7193 * ratio_spacing_height(icir) - 0.061) + &
664	!!$ (-1.273 * ratio_spacing_height(icir) + 0.9701) * ((1.1127 * &
665	!!$ ratio_spacing_height(icir) + 0.0311)*0.0) max_gap_factor(icir)
666	!!$
667	!!$ ! Calculating term_a and term_b for the new gust method (there are two methods):
668	!!$ IF (-2.1 * ratio_spacing_height(icir) + 0.91 < zero) THEN
669	!!$ term_a(icir) = zero
670	!!$ ELSE
671	!!$ term_a(icir) = -2.1 * ratio_spacing_height(icir) + 0.91
672	!!$ END IF
673	!!$ term_b(icir) = 1.0611 * LOG(ratio_spacing_height(icir)) + 4.2
674	!!$
675	!!$ ! Old gust factors:
676	!!$ old_gust_closer(icir) = max_bm_closer(icir) / mean_bm_closer(icir)
677	!!$ old_gust_further(icir) = max_bm_further(icir) / mean_bm_further(icir)
678	!!$ old_gust_edge(icir) = max_bm_edge(icir) / mean_bm_edge(icir)
679	!!$ old_gust_edge_gap(icir) = max_bm_edge_gap(icir) / mean_bm_edge_gap(icir)
680	!!$
681	!!$ ! New gust factors:
682	!!$ new_gust_closer(icir) = term_a(icir) * tree_heights_from_edge(icir) + &
683	!!$ term_b(icir)
684	!!$ new_gust_further(icir) = term_a(icir) * 9.0 + term_b(icir)
685	!!$ new_gust_edge(icir) = term_a(icir) * zero + term_b(icir)
686	!!$
687	!!$ ! Calculating the new gustiness at the edge accounting for the gap. This is done using the
688	!!$ ! argument that the ratios between old and new gustiness values stay the same:
689	!!$ IF ((old_gust_closer(icir) - old_gust_edge(icir)) > zero) THEN
690	!!$ new_gust_edge_gap_closer(icir) = (old_gust_edge_gap(icir) - &
691	!!$ old_gust_edge(icir)) * (new_gust_closer(icir) - &
692	!!$ new_gust_edge(icir)) / (old_gust_closer(icir) - &
693	!!$ old_gust_edge(icir)) + new_gust_edge(icir)
694	!!$ ELSE
695	!!$ new_gust_edge_gap_closer(icir) = new_gust_edge(icir)
696	!!$ END IF
697	!!$
698	!!$ IF ((old_gust_further(icir) - old_gust_edge(icir)) > zero) THEN
699	!!$ new_gust_edge_gap_further(icir) = (old_gust_edge_gap(icir) - &
700	!!$ old_gust_edge(icir)) * (new_gust_further(icir) - &
701	!!$ new_gust_edge(icir)) / (old_gust_further(icir) - &
702	!!$ old_gust_edge(icir)) + new_gust_edge(icir)
703	!!$ ELSE
704	!!$ new_gust_edge_gap_further(icir) = new_gust_edge(icir)
705	!!$ END IF
706	!!$
707	!!$ ! Calculating gustiness (g):
708	!!$ g_closer(icir) = ((mean_height(icir) * tree_heights_from_edge(icir) / &
709	!!$ (mean_height(icir) * 9.0)) * (new_gust_closer(icir) - new_gust_edge_gap_closer(icir))) +&
710	!!$ new_gust_edge_gap_closer(icir)
711	!!$
712	!!$ g_further(icir) = ((mean_height(icir) * 9.0 / (mean_height(icir) * 9.0)) * &
713	!!$ (new_gust_further(icir) - new_gust_edge_gap_further(icir))) +&
714	!!$ new_gust_edge_gap_further(icir)
715	!!$
716	!!$ ! The following section is used to calculate the edge_factor. This section calculates
717	!!$ ! how the mean loading on a tree changes as a function of tree height, spacing, distance
718	!!$ ! from edge and size of any upwind gap. An upwind gap greater than 10 tree heights
719	!!$ ! is assumed to be an infinite gap. The output (edge_factor) is used to modify the
720	!!$ ! calculated wind loading on the tree, which assumes a steady wind and a position well
721	!!$ ! inside the forest.
722	!!$ mean_bm_gust_factor_closer(icir) = (0.68 * ratio_spacing_height(icir) - &
723	!!$ 0.0385) + (-0.68 * ratio_spacing_height(icir) + 0.4785) * ((1.7239 * &
724	!!$ ratio_spacing_height(icir) + 0.0316)*tree_heights_from_edge(icir)) &
725	!!$ mean_gap_factor(icir)
726	!!$
727	!!$ mean_bm_gust_factor_further(icir) = (0.68 * ratio_spacing_height(icir) - &
728	!!$ 0.0385) + (-0.68 * ratio_spacing_height(icir) + 0.4785) * ((1.7239 * &
729	!!$ ratio_spacing_height(icir) + 0.0316)*9.0) mean_gap_factor(icir)
730	!!$
731	!!$ edge_factor(icir) = mean_bm_gust_factor_closer(icir) / &
732	!!$ mean_bm_gust_factor_further(icir)
733	!!$
734	!!$ ! Calculating gust_factor acounting for gustiness and edge_factor:
735	!!$ gust_factor_closer(icir) = g_closer(icir) * edge_factor(icir)
736	!!$ gust_factor_further(icir) = g_further(icir) * edge_factor(icir)
737	!!$
738	!!$
739	!!$
740	!!$ !! 5. Calculating critical moments (tree mechanics section)
741	!!$
742	!!$ ! overturning_moment_multiplier: this comes form the species parameters file and is related to the
743	!!$ ! depth and type of the soil.
744	!!$ ! Rooting depth has three categories: Shallow (less then 0.8 m deep); Deep (deeper than 0.8 m);
745	!!$ ! Average (average values to be used when rooting depth is unknown).
746	!!$ ! Soil type has four categories: Free-draining mineral soils; Gleyed mineral soils; Peaty mineral
747	!!$ ! soils; Deep peats
748	!!$ ! As these are generic soil types, a modified soil map of Europe would be needed for the optimal use
749	!!$ ! of WINDFALL.
750	!!$
751	!!$ IF (soil_type(ipts) == 'free_draining' .AND. rooting_depth(ipts,ivm) == 'shallow' &
752	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
753	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_shallow
754	!!$
755	!!$ ELSE IF (soil_type(ipts) == 'free_draining' .AND. rooting_depth(ipts,ivm) == 'shallow' &
756	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
757	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_shallow_leafless
758	!!$
759	!!$ ELSE IF (soil_type(ipts) == 'free_draining' .AND. rooting_depth(ipts,ivm) == 'deep' &
760	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
761	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_deep
762	!!$
763	!!$ ELSE IF (soil_type(ipts) == 'free_draining' .AND. rooting_depth(ipts,ivm) == 'deep' &
764	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
765	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_deep_leafless
766	!!$
767	!!$ ELSE IF (soil_type(ipts) == 'free_draining' .AND. rooting_depth(ipts,ivm) == 'average' &
768	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
769	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_average
770	!!$
771	!!$ ELSE IF (soil_type(ipts) == 'free_draining' .AND. rooting_depth(ipts,ivm) == 'average' &
772	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
773	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_free_draining_average_leafless
774	!!$
775	!!$ ELSE IF (soil_type(ipts) == 'gleyed' .AND. rooting_depth(ipts,ivm) == 'shallow' &
776	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
777	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_shallow
778	!!$
779	!!$ ELSE IF (soil_type(ipts) == 'gleyed' .AND. rooting_depth(ipts,ivm) == 'shallow' &
780	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
781	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_shallow_leafless
782	!!$
783	!!$ ELSE IF (soil_type(ipts) == 'gleyed' .AND. rooting_depth(ipts,ivm) == 'deep' &
784	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
785	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_deep
786	!!$
787	!!$ ELSE IF (soil_type(ipts) == 'gleyed' .AND. rooting_depth(ipts,ivm) == 'deep' &
788	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
789	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_deep_leafless
790	!!$
791	!!$ ELSE IF (soil_type(ipts) == 'gleyed' .AND. rooting_depth(ipts,ivm) == 'average' &
792	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
793	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_average
794	!!$
795	!!$ ELSE IF (soil_type(ipts) == 'gleyed' .AND. rooting_depth(ipts,ivm) == 'average' &
796	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
797	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_gleyed_average_leafless
798	!!$
799	!!$ ELSE IF (soil_type(ipts) == 'peaty' .AND. rooting_depth(ipts,ivm) == 'shallow' &
800	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
801	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peaty_shallow
802	!!$
803	!!$ ELSE IF (soil_type(ipts) == 'peaty' .AND. rooting_depth(ipts,ivm) == 'shallow' &
804	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
805	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peaty_shallow_leafless
806	!!$
807	!!$ ELSE IF (soil_type(ipts) == 'peaty' .AND. rooting_depth(ipts,ivm) == 'deep' &
808	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
809	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peaty_deep
810	!!$
811	!!$ ELSE IF (soil_type(ipts) == 'peaty' .AND. rooting_depth(ipts,ivm) == 'deep' &
812	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
813	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peaty_deep_leafless
814	!!$
815	!!$ ELSE IF (soil_type(ipts) == 'peaty' .AND. rooting_depth(ipts,ivm) == 'average' &
816	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
817	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peaty_average
818	!!$
819	!!$ ELSE IF (soil_type(ipts) == 'peaty' .AND. rooting_depth(ipts,ivm) == 'average' &
820	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
821	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peaty_average_leafless
822	!!$
823	!!$ ELSE IF (soil_type(ipts) == 'peaty' .AND. rooting_depth(ipts,ivm) == 'shallow' &
824	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
825	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peat_shallow
826	!!$
827	!!$ ELSE IF (soil_type(ipts) == 'peat' .AND. rooting_depth(ipts,ivm) == 'shallow' &
828	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
829	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peat_shallow_leafless
830	!!$
831	!!$ ELSE IF (soil_type(ipts) == 'peat' .AND. rooting_depth(ipts,ivm) == 'deep' &
832	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
833	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peat_deep
834	!!$
835	!!$ ELSE IF (soil_type(ipts) == 'peat' .AND. rooting_depth(ipts,ivm) == 'deep' &
836	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
837	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peat_deep_leafless
838	!!$
839	!!$ ELSE IF (soil_type(ipts) == 'peat' .AND. rooting_depth(ipts,ivm) == 'average' &
840	!!$ .AND. senescence(ipts,ivm) .EQ. .FALSE.) THEN
841	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peat_average
842	!!$
843	!!$ ELSE IF (soil_type(ipts) == 'peat' .AND. rooting_depth(ipts,ivm) == 'average' &
844	!!$ .AND. senescence(ipts,ivm) .EQ. .TRUE.) THEN
845	!!$ overturning_moment_multiplier(ipts,ivm) = overturning_peat_average_leafless
846	!!$
847	!!$ ELSE
848	!!$ overturning_moment_multiplier(ipts,ivm) = 125.0
849	!!$ ! average value if parameters are missing
850	!!$ END IF
851	!!$
852	!!$ ! Maximum moment calculated for tree overturning.
853	!!$ max_overturning_moment(ipts,ivm,icir) = overturning_moment_multiplier(ipts,ivm) * &
854	!!$ stem_mass(ipts,ivm,icir)
855	!!$
856	!!$ ! Maximum moment calculated for stem breakage.
857	!!$ max_breaking_moment(ipts,ivm,icir) = modulus_rupture * f_knot * pi * &
858	!!$ (section_diam(ipts,ivm,icir,1)**3.0) / 32.0
859	!!$
860	!!$
861	!!$ ! This is where the gustiness of wind is "filtered out"
862	!!$ ! so that hourly mean wind speeds can be used in the simulations.
863	!!$ ! In the end, moments and hence, critical wind speeds are calculated
864	!!$ ! for both areas (closer & further) with different gustiness.
865	!!$
866	!!$ overturning_moment_closer(ipts,ivm,icir) = max_overturning_moment(ipts,ivm,icir) / &
867	!!$ gust_factor_closer(ipts,ivm,icir)
868	!!$
869	!!$ overturning_moment_further(ipts,ivm,icir) = max_overturning_moment(ipts,ivm,icir) / &
870	!!$ gust_factor_further(ipts,ivm,icir)
871	!!$
872	!!$ breaking_moment_closer(ipts,ivm,icir) = max_breaking_moment(ipts,ivm,icir) / &
873	!!$ gust_factor_closer(ipts,ivm,icir)
874	!!$
875	!!$ breaking_moment_further(ipts,ivm,icir) = max_breaking_moment(ipts,ivm,icir) / &
876	!!$ gust_factor_further(ipts,ivm,icir)
877
878	!!$
879	!!$
880	!!$ !! 7. Calculating critical wind speed and damage rate
881	!!$ !! 7.1 Calculating critical wind speed in the forest and 10 m above the zero-plane displacement
882	!!$ ! for the edge area
883	!!$
884	!!$ ! Critical wind speed for stem breakage in the edge area (m/s)
885	!!$ call calculate_wind_process(breaking_moment_closer(ipts,ivm,icir),cws_break_closer(ipts,ivm,icir))
886	!!$ ! Critical wind speed for stem breakage 10 m above the zero-plane displacement in the edge area
887	!!$ ! (m/s)
888	!!$ call elevate(cws_break_closer(ipts,ivm,icir), z0, d, cws_break_closer_10(ipts,ivm,icir))
889	!!$
890	!!$ ! Critical wind speed for tree overturning in the edge area (m/s)
891	!!$ call calculate_wind_process(overturning_moment_closer(ipts,ivm,icir),cws_overturn_closer(ipts,ivm,icir))
892	!!$ ! Critical wind speed for tree overturning 10 m above the zero-plane displacement in the edge area
893	!!$ ! (m/s)
894	!!$ call elevate(cws_overturn_closer(ipts,ivm,icir), z0, d, cws_overturn_closer_10)
895	!!$
896	!!$
897	!!$ ! Calculating damage severity in case any of the critical wind speeds is reached in the area
898	!!$ ! around gaps (forest edge area):
899	!!$
900	!!$ ! Until a better method is found, we kill all the trees in the given dbh-class of the PFT
901	!!$ ! in case of damage.
902	!!$
903	!!$ ! Whichever of the two damage types (overturning, stem breakage) has the lowest critical wind speed,
904	!!$ ! it will be the occurring one.
905	!!$ IF (cws_overturn_closer_10(ipts,ivm,icir) < cws_break_closer_10(ipts,ivm,icir)) THEN
906	!!$ cws_final_closer(ipts,ivm,icir) = cws_overturn_closer_10(ipts,ivm,icir)
907	!!$ ELSE
908	!!$ cws_final_closer(ipts,ivm,icir) = cws_break_closer_10(ipts,ivm,icir)
909	!!$ END IF
910	!!$
911	!!$ IF (wind_speed_actual(ipts) >= cws_final_closer(ipts,ivm,icir)) THEN
912	!!$ IF (cws_overturn_closer_10(ipts,ivm,icir) < cws_break_closer_10(ipts,ivm,icir)) THEN
913	!!$ wind_damage_type_closer(ipts,ivm,icir) = overturning
914	!!$ ELSE
915	!!$ wind_damage_type_closer(ipts,ivm,icir) = breakage
916	!!$ END IF
917	!!$ wind_damage_rate_closer(ipts,ivm,icir) = 1.0
918	!!$ ELSE
919	!!$ wind_damage_rate_closer(ipts,ivm,icir) = zero
920	!!$ ! Damage rates should be dependent on the actual difference between cws and actual wind speed.
921	!!$ ! However, as literature is lacking, we use 100% mortality rates at the moment (as in GALES).
922	!!$ END IF
923	!!$
924	!!$
925	!!$ !! 7.2 Calculating critical wind speed in the forest and 10 m above the zero-plane displacement
926	!!$ ! for the inner forest area
927	!!$
928	!!$ ! These are only calculated if there was a damage rate higher than zero in the area around the gaps
929	!!$ ! (forest edge area).
930	!!$
931	!!$ IF (wind_damage_rate_closer(ipts,ivm,icir) > zero) THEN
932	!!$
933	!!$ ! Whichever of the two damage types (overturning, stem breakage) has the lowest critical wind
934	!!$ ! speed, it will be the occurring one.
935	!!$
936	!!$ ! Critical wind speed for stem breakage in the inner forest area (m/s)
937	!!$ call calculate_wind_process(breaking_moment_further(ipts,ivm,icir),cws_break_further(ipts,ivm,icir))
938	!!$ ! Critical wind speed for stem breakage 10 m above the zero-plane displacement in the inner
939	!!$ ! forest area (m/s)
940	!!$ call elevate(cws_break_further(ipts,ivm,icir), surf_rough(ipts,ivm,icir), &
941	!!$ zp_displ(ipts,ivm,icir), cws_break_further_10(ipts,ivm,icir))
942	!!$
943	!!$ ! Critical wind speed for tree overturning in the inner forest area (m/s)
944	!!$ call calculate_wind_process(overturning_moment_further(ipts,ivm,icir),&
945	!!$ cws_overturn_further(ipts,ivm,icir))
946	!!$ ! Critical wind speed for tree overturning 10 m above the zero-plane displacement in the
947	!!$ ! inner forest area (m/s)
948	!!$ call elevate(cws_overturn_further(ipts,ivm,icir), surf_rough(ipts,ivm,icir), &
949	!!$ zp_displ(ipts,ivm,icir), cws_overturn_further_10(ipts,ivm,icir))
950	!!$
951	!!$ ! Whichever of the two damage types (overturning, stem breakage) has the lowest critical wind
952	!!$ ! speed, it will be the occurring one.
953	!!$ IF (cws_overturn_further_10(ipts,ivm,icir) < cws_break_further_10(ipts,ivm,icir)) THEN
954	!!$ cws_final_further(ipts,ivm,icir) = cws_overturn_further_10(ipts,ivm,icir)
955	!!$ ELSE
956	!!$ cws_final_further(ipts,ivm,icir) = cws_break_further_10(ipts,ivm,icir)
957	!!$ END IF
958	!!$
959	!!$ IF (wind_speed_actual(ipts) >= cws_final_further(ipts,ivm,icir)) THEN
960	!!$ IF (cws_overturn_further_10(ipts,ivm,icir) < cws_break_closer_10(ipts,ivm,icir)) THEN
961	!!$ wind_damage_type_further(ipts,ivm,icir) = overturning
962	!!$ ELSE
963	!!$ wind_damage_type_further(ipts,ivm,icir) = breakage
964	!!$ END IF
965	!!$ wind_damage_rate_further(ipts,ivm,icir) = 1.0
966	!!$ ELSE
967	!!$ wind_damage_rate_further(ipts,ivm,icir) = zero
968	!!$ END IF
969	!!$
970	!!$ ELSE
971	!!$ wind_damage_rate_further(ipts,ivm,icir) = zero
972	!!$ END IF
973	!!$
974	!!$
975	!!$
976	!!$ ! NOTE: this is the point where damage should be taken into account in other
977	!!$ ! modules of ORCHIDEE. Output variables should be calculated here.
978	!!$
979	!!$
980	!!$
981	!!$ END DO ! loop ncirc
982
983	END DO ! loop nvm
984
985	END DO ! loop npts
986
987
988
989	!!$ !! 9. Writing out history variables
990	!!$
991	!!$ ! NOTE: THIS IS ADOPTED FROM stomate_fire AND HASN'T CHANGED YET.
992	!!$ firefrac(:,:) = firefrac(:,:) / dt
993	!!$
994	!!$ CALL xios_orchidee_send_field("FIREFRAC",firefrac(:,:))
995	!!$ CALL xios_orchidee_send_field("FIREDEATH",firedeath(:,:))
996	!!$
997	!!$ CALL histwrite_p (hist_id_stomate, 'FIREFRAC', itime, &
998	!!$ firefrac(:,:), npts*nvm, horipft_index)
999	!!$ CALL histwrite_p (hist_id_stomate, 'FIREDEATH', itime, &
1000	!!$ firedeath(:,:), npts*nvm, horipft_index)
1001	!!$
1002	!!$ IF (bavard.GE.4) WRITE(numout,*) 'Leaving fire'
1003	!!$ ELSE
1004	!!$ ! This is for bavard.
1005	!!$ END IF
1006
1007
1008
1009	END SUBROUTINE wind_damage
1010
1011	!!$
1012	!!$
1013	!!$!! ===========================================================================================
1014	!!$!! SUBROUTINE : streamlining
1015	!!$!!
1016	!!$!>\BRIEF : This subroutine calculates the streamlining of tree crowns. Streamlining is
1017	!!$!! the change of shape of the crowns due to wind.
1018	!!$!!
1019	!!$!! DESCRIPTION : None
1020	!!$!!
1021	!!$!! RECENT CHANGE(S): None
1022	!!$!!
1023	!!$!! RETURN VALUE : Surface roughness and zero-plane displacement among other variables that
1024	!!$!! are used in the other subroutines
1025	!!$!! REFERENCES :
1026	!!$!!
1027	!!$!! FLOWCHART : None
1028	!!$!! \streamlining_n
1029	!!$!_ ============================================================================================
1030	!!$
1031	!!$ SUBROUTINE streamlining(a_guess_wind_speed, porosity_sub, lambda_sub, lambda_capital_sub, &
1032	!!$ gamma_solved_sub, zp_displ_sub, psih_sub, surf_rough_sub)
1033	!!$
1034	!!$ IMPLICIT NONE
1035	!!$
1036	!!$!! 0. Variable and parameter declaration
1037	!!$
1038	!!$ !! 0.1 Input variables
1039	!!$
1040	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: a_guess_wind_speed
1041	!!$
1042	!!$ !! 0.2 Output variables
1043	!!$
1044	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: porosity_sub
1045	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: lambda_sub
1046	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: lambda_capital_sub
1047	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: gamma_solved_sub
1048	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: zp_displ_sub
1049	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: psih_sub
1050	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: surf_rough_sub
1051	!!$
1052	!!$ !! 0.3 Modified variables
1053	!!$
1054	!!$ !! 0.4 Local variables
1055	!!$
1056	!!$!_ =============================================================================================
1057	!!$
1058	!!$!! 1. Subroutine
1059	!!$
1060	!!$ IF (a_guess_wind_speed(ipts,ivm,icir) > 25.0) THEN
1061	!!$ ! Above 25 m/s, tree crowns cannot streamline any further.
1062	!!$ porosity_sub(ipts,ivm,icir) = streamlining_c * 25.0**(-streamlining_n)
1063	!!$ ELSE
1064	!!$ porosity_sub(ipts,ivm,icir) = streamlining_c * &
1065	!!$ a_guess_wind_speed(ipts,ivm,icir)**(-streamlining_n)
1066	!!$ END IF
1067	!!$
1068	!!$ lambda_sub(ipts,ivm,icir) = (canopy_breadth(ipts,ivm,icir)/2.0) * canopy_depth(ipts,ivm,icir) * &
1069	!!$ porosity_sub(ipts,ivm,icir) / current_spacing(ipts,ivm,icir)**2.0
1070	!!$ lambda_capital_sub(ipts,ivm,icir) = 2.0 * lambda_sub(ipts,ivm,icir)
1071	!!$
1072	!!$ IF (lambda_capital_sub(ipts,ivm,icir) > 0.6) THEN
1073	!!$ ! Gamma is needed for the calculation of the roughness length.
1074	!!$ gamma_solved_sub(ipts,ivm,icir) = 1.0 / ((c_surface + c_drag * 0.3)**0.5)
1075	!!$ ELSE
1076	!!$ gamma_solved_sub(ipts,ivm,icir) = 1.0 / ((c_surface + c_drag * &
1077	!!$ lambda_capital_sub(ipts,ivm,icir) / 2.0)**0.5)
1078	!!$ END IF
1079	!!$
1080	!!$ ! Calculation of the zero-plane displacement:
1081	!!$ zp_displ_sub(ipts,ivm,icir) = mean_height(ipts,ivm,icir) * (1.0 - ((1.0 - EXP(-(c_displacement * &
1082	!!$ lambda_capital_sub(ipts,ivm,icir))**0.5)) / &
1083	!!$ (c_displacement * lambda_capital_sub(ipts,ivm,icir))**0.5))
1084	!!$
1085	!!$ ! Calculation of the aerodynamic roughness length:
1086	!!$
1087	!!$ ! psih = "Profile influence function" at h (height of the roughness elements).
1088	!!$ psih_sub(ipts,ivm,icir) = LOG(c_roughness) -1.0 + c_roughness**(-1.0)
1089	!!$ surf_rough_sub(ipts,ivm,icir) = (mean_height(ipts,ivm,icir) - zp_displ_sub(ipts,ivm,icir)) * &
1090	!!$ EXP(psih_sub(ipts,ivm,icir) - ct_karman * gamma_solved_sub(ipts,ivm,icir))
1091	!!$
1092	!!$ END SUBROUTINE streamlining
1093	!!$
1094	!!$
1095	!!$
1096	!!$!! ============================================================================================
1097	!!$!! SUBROUTINE : calculate_force
1098	!!$!!
1099	!!$!>\BRIEF : This subroutine calculates the force affecting the tree.
1100	!!$!!
1101	!!$!! DESCRIPTION : None
1102	!!$!!
1103	!!$!! RECENT CHANGE(S): None
1104	!!$!!
1105	!!$!! RETURN VALUE :
1106	!!$!!
1107	!!$!! REFERENCES :
1108	!!$!!
1109	!!$!! FLOWCHART : None
1110	!!$!! \streamlining_n
1111	!!$!_ ============================================================================================
1112	!!$
1113	!!$ SUBROUTINE calculate_force(a_guess_wind_speed, a_gamma_solved, a_zp_displ, force_on_tree_sub, &
1114	!!$ height_of_force_sub)
1115	!!$
1116	!!$ IMPLICIT NONE
1117	!!$
1118	!!$!! 0. Variable and parameter declaration
1119	!!$
1120	!!$ !! 0.1 Input variables
1121	!!$
1122	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: a_guess_wind_speed
1123	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: a_gamma_solved
1124	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: a_zp_displ
1125	!!$
1126	!!$ !! 0.2 Output variables
1127	!!$
1128	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: force_on_tree_sub
1129	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: height_of_force_sub
1130	!!$
1131	!!$ !! 0.3 Modified variables
1132	!!$
1133	!!$ !! 0.4 Local variables
1134	!!$
1135	!!$!_ =============================================================================================
1136	!!$
1137	!!$!! 1. Subroutine
1138	!!$
1139	!!$ force_on_tree_sub(ipts,ivm,icir) = air_density * (a_guess_wind_speed(ipts,ivm,icir) * &
1140	!!$ current_spacing(ipts,ivm,icir) / a_gamma_solved(ipts,ivm,icir))**2.0
1141	!!$
1142	!!$ height_of_force_sub(ipts,ivm,icir) = a_zp_displ(ipts,ivm,icir)
1143	!!$
1144	!!$ ! The maximum of the zero-plane displacement height is locked at 0.8 tree height.
1145	!!$ IF (height_of_force_sub(ipts,ivm,icir) > 0.8 * mean_height(ipts,ivm,icir)) THEN
1146	!!$ height_of_force_sub(ipts,ivm,icir) = 0.8 * mean_height(ipts,ivm,icir)
1147	!!$ END IF
1148	!!$
1149	!!$ END subroutine calculate_force
1150	!!$
1151	!!$
1152	!!$!! =============================================================================================
1153	!!$!! SUBROUTINE : calculate_bending_moment
1154	!!$!!
1155	!!$!>\BRIEF : Subroutine to calculate the bending moment
1156	!!$!! (to compare with critical moments in the iterative calculation of the critical
1157	!!$!! wind speeds)
1158	!!$!!
1159	!!$!! DESCRIPTION : None
1160	!!$!!
1161	!!$!! RECENT CHANGE(S): None
1162	!!$!!
1163	!!$!! RETURN VALUE :
1164	!!$!!
1165	!!$!! REFERENCES :
1166	!!$!!
1167	!!$!! FLOWCHART : None
1168	!!$!! \streamlining_n
1169	!!$!_ =============================================================================================
1170	!!$
1171	!!$ SUBROUTINE calculate_bending_moment (a_force_on_tree, a_height_of_force, bending_moment_sub)
1172	!!$
1173	!!$ IMPLICIT NONE
1174	!!$
1175	!!$!! 0. Variable and parameter declaration
1176	!!$
1177	!!$ !! 0.1 Input variables
1178	!!$
1179	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: a_force_on_tree
1180	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: a_height_of_force
1181	!!$
1182	!!$ !! 0.2 Output variables
1183	!!$
1184	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: bending_moment_sub
1185	!!$
1186	!!$ !! 0.3 Modified variables
1187	!!$
1188	!!$ !! 0.4 Local variables
1189	!!$
1190	!!$!_ =============================================================================================
1191	!!$
1192	!!$!! 1. Subroutine
1193	!!$
1194	!!$ bending_moment_sub(ipts,ivm,icir) = a_force_on_tree(ipts,ivm,icir) * &
1195	!!$ a_height_of_force(ipts,ivm,icir) * f_crown_mass
1196	!!$
1197	!!$ END SUBROUTINE calculate_bending_moment
1198	!!$
1199	!!$
1200	!!$
1201	!!$!! =============================================================================================
1202	!!$!! SUBROUTINE : calculate_cws
1203	!!$!!
1204	!!$!>\BRIEF : Subroutine to calculate critical wind speed
1205	!!$!! (to compare with wind speeds used for calculating streamlining in the
1206	!!$!! iterative calculation of the final critical wind speed)
1207	!!$!!
1208	!!$!! DESCRIPTION : None
1209	!!$!!
1210	!!$!! RECENT CHANGE(S): None
1211	!!$!!
1212	!!$!! RETURN VALUE :
1213	!!$!!
1214	!!$!! REFERENCES :
1215	!!$!!
1216	!!$!! FLOWCHART : None
1217	!!$!! \streamlining_n
1218	!!$!_ =============================================================================================
1219	!!$
1220	!!$ SUBROUTINE calculate_cws(critical_moment_sub, zp_displ_sub, gamma_solved_sub, &
1221	!!$ current_spacing_sub, crit_wind_speed_sub)
1222	!!$
1223	!!$ IMPLICIT NONE
1224	!!$
1225	!!$!! 0. Variable and parameter declaration
1226	!!$
1227	!!$ !! 0.1 Input variables
1228	!!$
1229	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: critical_moment_sub
1230	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: zp_displ_sub
1231	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: gamma_solved_sub
1232	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: current_spacing_sub
1233	!!$ ! f_crown_mass and rho are parameters.
1234	!!$
1235	!!$ !! 0.2 Output variables
1236	!!$
1237	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: crit_wind_speed_sub
1238	!!$
1239	!!$ !! 0.3 Modified variables
1240	!!$
1241	!!$ !! 0.4 Local variables
1242	!!$
1243	!!$!_ =============================================================================================
1244	!!$
1245	!!$!! 1. Subroutine
1246	!!$
1247	!!$ crit_wind_speed_sub(ipts,ivm,icir) = (gamma_solved_sub(ipts,ivm,icir) / &
1248	!!$ current_spacing_sub(ipts,ivm,icir)) * (critical_moment_sub(ipts,ivm,icir) / &
1249	!!$ (f_crown_mass * air_density * zp_displ_sub(ipts,ivm,icir)))**0.5
1250	!!$
1251	!!$ END SUBROUTINE calculate_cws
1252	!!$
1253	!!$
1254	!!$
1255	!!$!! =============================================================================================
1256	!!$!! SUBROUTINE : calculate_wind_process
1257	!!$!!
1258	!!$!>\BRIEF : This subroutine calculates the breaking and overturning critical wind
1259	!!$!! speeds
1260	!!$!!
1261	!!$!! DESCRIPTION : None
1262	!!$!!
1263	!!$!! RECENT CHANGE(S): None
1264	!!$!!
1265	!!$!! RETURN VALUE :
1266	!!$!!
1267	!!$!! REFERENCES :
1268	!!$!!
1269	!!$!! FLOWCHART : None
1270	!!$!! \streamlining_n
1271	!!$!_ =============================================================================================
1272	!!$
1273	!!$ SUBROUTINE calculate_wind_process(critical_moment, guess_wind_speed_sub)
1274	!!$
1275	!!$ IMPLICIT NONE
1276	!!$
1277	!!$!! 0. Variable and parameter declaration
1278	!!$
1279	!!$ !! 0.1 Input variables
1280	!!$
1281	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: critical_moment
1282	!!$
1283	!!$ !! 0.2 Output variables
1284	!!$
1285	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: guess_wind_speed_sub !! A guess for the
1286	!!$ !! critical wind speed
1287	!!$
1288	!!$ !! 0.3 Modified variables
1289	!!$
1290	!!$ !! 0.4 Local variables
1291	!!$
1292	!!$ REAL, DIMENSION(npts,nvm,ncirc) :: delta !! Difference between the guessed
1293	!!$ !! and the calculated wind speed
1294	!!$ REAL, DIMENSION(npts,nvm,ncirc) :: wind_precision !! How close we should
1295	!!$ !! get with the iterations
1296	!!$
1297	!!$!_ =============================================================================================
1298	!!$
1299	!!$!! 1. Subroutine
1300	!!$
1301	!!$! This bit belongs to the original form of doing the iterations. For more info, see the notes below.
1302	!!$! guess_wind_speed_sub = 64.0
1303	!!$! delta = guess_wind_speed_sub / 2.0
1304	!!$
1305	!!$ guess_wind_speed_sub(ipts,ivm,icir) = 25.0
1306	!!$
1307	!!$ wind_precision(ipts,ivm,icir) = 0.1
1308	!!$
1309	!!$ DO WHILE (delta(ipts,ivm,icir) > wind_precision(ipts,ivm,icir))
1310	!!$ CALL streamlining(guess_wind_speed_sub(ipts,ivm,icir), porosity(ipts,ivm,icir), &
1311	!!$ lambda(ipts,ivm,icir), lambda_capital(ipts,ivm,icir), gamma_solved(ipts,ivm,icir), &
1312	!!$ zp_displ(ipts,ivm,icir), psih(ipts,ivm,icir), surf_rough(ipts,ivm,icir))
1313	!!$
1314	!!$! The following lines belong to the original way of calculating the final critical wind speed
1315	!!$! with multiple iterations, as streamlining depends on wind speed and critical wind speed
1316	!!$! depends on streamlining. Originally, iterations started with a wind speed value of 64 m/s and
1317	!!$! then if bending moment was below the critical moment, then a wind speed value of 32 m/s was
1318	!!$! used. If bending moment was above critical moment afterwards, then a wind speed value of 48
1319	!!$! m/s was used, and so on, until the difference between the critical moment and bending moment
1320	!!$! was less than the pre-defined precision.
1321	!!$! The new version uses a shorter iterating process instead.
1322	!!$! Now the critical wind speed is calculated instead of the moments, and the iterating process
1323	!!$! starts from 25 m/s (no further streamlining above that value) and uses the calculated critical
1324	!!$! wind speed as the wind speed used for calculating streamlining after that first value.
1325	!!$! Then, with the new streamlining values a new critical wind speed is calculated, and so on,
1326	!!$! until the difference between the two wind speeds (guess_wind_speed_sub and crit_wind_speed)
1327	!!$! is less than the pre-defined precision.
1328	!!$
1329	!!$! call calculate_force(guess_wind_speed_sub, gammaSolved, d, force_on_tree, height_of_force)
1330	!!$! gammaSolved and d come from call streamlining.
1331	!!$! call calculate_bending_moment(force_on_tree, height_of_force, bending_moment)
1332	!!$! force_on_tree and height_of_force come from call calculate_force.
1333	!!$! IF (bending_moment > critical_moment) THEN
1334	!!$! guess_wind_speed_sub = guess_wind_speed_sub - delta
1335	!!$! ELSE
1336	!!$! guess_wind_speed_sub = guess_wind_speed_sub + delta
1337	!!$! END IF
1338	!!$! delta = delta / 2.0
1339	!!$! END DO
1340	!!$
1341	!!$ ! The input variable of the subroutine is passed towards the call of another subroutine:
1342	!!$ critical_moment_cws(ipts,ivm,icir) = critical_moment(ipts,ivm,icir)
1343	!!$
1344	!!$ CALL calculate_cws(critical_moment_cws(ipts,ivm,icir), zp_displ(ipts,ivm,icir), &
1345	!!$ gamma_solved(ipts,ivm,icir), current_spacing(ipts,ivm,icir), &
1346	!!$ crit_wind_speed(ipts,ivm,icir))
1347	!!$
1348	!!$ delta(ipts,ivm,icir) = guess_wind_speed_sub(ipts,ivm,icir) - &
1349	!!$ crit_wind_speed(ipts,ivm,icir)
1350	!!$
1351	!!$ guess_wind_speed_sub(ipts,ivm,icir) = crit_wind_speed(ipts,ivm,icir)
1352	!!$ END DO
1353	!!$
1354	!!$ END SUBROUTINE calculate_wind_process
1355	!!$
1356	!!$
1357	!!$
1358	!!$!! =============================================================================================
1359	!!$!! SUBROUTINE : elevate
1360	!!$!!
1361	!!$!>\BRIEF : This subroutine converts the critical wind speeds calculated above to the
1362	!!$!! wind speed 10 meters above the zero-plane displacement
1363	!!$!!
1364	!!$!! DESCRIPTION : None
1365	!!$!!
1366	!!$!! RECENT CHANGE(S): None
1367	!!$!!
1368	!!$!! RETURN VALUE :
1369	!!$!!
1370	!!$!! REFERENCES :
1371	!!$!!
1372	!!$!! FLOWCHART : None
1373	!!$!! \streamlining_n
1374	!!$!_ =============================================================================================
1375	!!$
1376	!!$ SUBROUTINE elevate(uh_speed, the_surf_rough, the_zp_displ, elevate_result_sub)
1377	!!$
1378	!!$ IMPLICIT NONE
1379	!!$
1380	!!$!! 0. Variable and parameter declaration
1381	!!$
1382	!!$ !! 0.1 Input variables
1383	!!$
1384	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: uh_speed
1385	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: the_surf_rough
1386	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(in) :: the_zp_displ
1387	!!$
1388	!!$ !! 0.2 Output variables
1389	!!$
1390	!!$ REAL, DIMENSION(npts,nvm,ncirc), INTENT(out) :: elevate_result_sub
1391	!!$
1392	!!$ !! 0.3 Modified variables
1393	!!$
1394	!!$ !! 0.4 Local variables
1395	!!$
1396	!!$!_ ==============================================================================================
1397	!!$
1398	!!$!! 1. Subroutine
1399	!!$
1400	!!$ elevate_result_sub(ipts,ivm,icir) = uh_speed(ipts,ivm,icir) * log(10 / &
1401	!!$ the_surf_rough(ipts,ivm,icir)) / log((mean_height(ipts,ivm,icir) - the_zp_displ(ipts,ivm,icir)) / &
1402	!!$ the_surf_rough(ipts,ivm,icir))
1403	!!$
1404	!!$ END SUBROUTINE elevate
1405
1406
1407
1408	END MODULE stomate_windfall

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

Download in other formats: