Context Navigation

stomate_laieff.f90 @ 7475

Last change on this file since 7475 was 2945, checked in by sebastiaan.luyssaert, 9 years ago

DEV: tested 1 year global. This code contains the latest version for anthropogenic tree species channges, several bug fixes to forest management as well as the code for the fully integrated multi-layer energy budget. This implies that the multi-layer energy budget makes use Pinty's albedo scheme, the rognostic canopy structure as well as a vertical profile for stomatal conductance. This is an intermediate version because species change code is not complete as some management changes have not been implemented yet. Further the multi-layer albedo code needs more work in terms of calculating average fluxes at the pixel rather than the PFT level

File size: 125.7 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_laieff
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Groups the subroutines that are used to calculate the effective
10	!! LAI. Right now this is only used in the albedo routines. I
11	!! originally tried to put these routines there, but that creates
12	!! some dependency problems that could not be resolved, since
13	!! the effective LAI is calculated in stomate.f90 along with
14	!! the LAI.
15	!!
16	!!\n DESCRIPTION : None
17	!!
18	!! RECENT CHANGE(S) : None
19	!!
20	!! REFERENCE(S) : None
21	!!
22	!! SVN :
23	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/stomate_laieff.f90 $
24	!! $Date: 2013-02-13 11:14:38 +0100 (Wed, 13 Feb 2013) $
25	!! $Revision: 1182 $
26	!! \n
27	!_ ================================================================================================================================
28
29	MODULE stomate_laieff
30
31	! Modules used:
32	USE constantes
33	USE pft_parameters
34	USE structures
35	USE ioipsl_para, ONLY : ipslerr_p
36	USE stomate_stand_structure
37	USE stomate_data
38	USE function_library, ONLY: wood_to_height_eff, wood_to_cn_eff, biomass_to_lai
39	USE matrix_resolution
40
41	IMPLICIT NONE
42
43	! Private & public routines
44
45	PUBLIC effective_lai, fitting_laieff, calculate_z_level_photo, &
46	find_lai_per_level, combine_lai_levels
47
48	INTEGER,SAVE :: ivm_save
49	INTEGER,SAVE :: ibin_save
50	INTEGER,SAVE :: ilevel_save
51
52	CONTAINS
53
54
55	!! ================================================================================================================================
56	!! SUBROUTINE : calculate_z_level_photo
57	!!
58	!>\BRIEF
59	!!
60	!! DESCRIPTION : This is routine to calculate the physical heights of
61	!! all the levels used to calculate photosynthesis. The bottom layer
62	!! will always have a height of zero, since it must be at the ground
63	!! level. I tried to make all the layers evenly spaced starting from
64	!! zero to the top of the canopy, but this does not work very well
65	!! for tall trees; all the vegetation ends up in the top level and
66	!! nothing in the lower levels. So now I try to find the bottom of
67	!! the canopy for trees as well, and start the second level from there.
68	!! The top level is always just below the top of the canopy (unless there
69	!! is only one level). If there are multiple levels used in the energy
70	!! budget, the spacing between the levels will not be uniform.
71	!!
72	!! If there is no leaf biomass, all the levels are set to zero since
73	!! they should not be used and it would be a waste of time to calculate.
74	!!
75	!! NOTE:
76	!!
77	!! RECENT CHANGE(S) : None
78	!!
79	!! MAIN OUTPUT VARIABLE(S): ::z_level_photo
80	!!
81	!! REFERENCE(S) :
82	!!
83	!! FLOWCHART : None
84	!! \n
85	!_ ================================================================================================================================
86
87	SUBROUTINE calculate_z_level_photo(npts, circ_class_biomass, &
88	circ_class_n, biomass, z_level_photo)
89
90	!! 0 Variable and parameter declaration
91
92	!! 0.1 Input variables
93
94	INTEGER(i_std),INTENT(IN) :: npts !! Domain size - number of pixels (unitless)
95	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(IN) &
96	:: circ_class_biomass !! Biomass of the componets of the model
97	!! tree within a circumference
98	!! class @tex $(gC ind^{-1})$ @endtex
99	REAL(r_std), DIMENSION(:,:,:), INTENT(IN) :: circ_class_n !! Number of trees within each circ class
100	!! class @tex $(m^{-2})$ @endtex
101	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! Stand level biomass @tex $(gC m^{-2})$ @endtex
102	!! 0.2 Output variables
103	REAL(r_std),DIMENSION(:,:,:),INTENT(OUT) :: z_level_photo !! the height of the bottom of each of our levels
104	!! @tex $(m)$ @endtex
105	!! note that level 1 is closest to the ground
106
107
108	!! 0.3 Modified variables
109
110	!! 0.4 Local variables
111	INTEGER :: ipts,ivm,il,il_photo,nl_photo,icir !! Indices (unitless)
112	REAL(r_std),DIMENSION(ncirc) :: heights !! Tree heights @tex $(m)$ @endtex
113	REAL(r_std),DIMENSION(ncirc) :: cn_dia !! Crown diameters @tex $(m)$ @endtex
114	REAL(r_std),DIMENSION(ncirc) :: cn_bottom !! Height of the bottom of the crown @tex $(m)$ @endtex
115	REAL(r_std) :: max_height !! Highest actual vegetation height @tex $(m)$ @endtex
116	REAL(r_std) :: min_height !! Lowest bottom of the canopy @tex $(m)$ @endtex
117	REAL(r_std) :: upper_level_boundary !! The height of the top of the current energy budget
118	!! level @tex $(m)$ @endtex
119	REAL(r_std) :: lower_level_boundary !! The height of the bottom of the current energy budget
120	!! level @tex $(m)$ @endtex
121	REAL(r_std) :: ideal_level_thickness !! The thickness of the photosynthesis levels in this
122	!! energy budget level if there are no other concerns.
123	!! @tex $(m)$ @endtex
124	REAL(r_std) :: block_thickness !! The thickness of a chunk of PS levels together.
125	!! @tex $(m)$ @endtex
126	REAL(r_std) :: canopy_thickness !! The thickness of the whole vegetation layer.
127	!! @tex $(m)$ @endtex
128	REAL(r_std) :: canopy_center !! The height of the center of the whole vegetation layer.
129	!! @tex $(m)$ @endtex
130	!_ ================================================================================================================================
131
132	IF (bavard.GE.2) WRITE(numout,*) 'Entering calculate_z_level_photo',nlevels,nlevels_tot
133
134	nl_photo=nlevels_tot/nlevels ! the number of photosynthesis levels per energy level
135
136	! If we are not using any additional levels for
137	! the photosynthesis, this is easy.
138	IF(nlevels_photo .EQ. 1) THEN
139	DO il=1,nlevels
140	z_level_photo(:,:,il)=z_level(il)
141	ENDDO
142	RETURN
143	ENDIF
144
145	IF(nlevels .GT. 1)THEN
146	WRITE(numout,*) 'WARNING: calculate_z_level_photo has not been tested yet'
147	WRITE(numout,*) 'WARNING: for the multilevel energy budget. Be sure to '
148	WRITE(numout,*) 'WARNING: do this and then remove this warning message.'
149	ENDIF
150
151	! This is just a quick test. The difference between the energy level heights
152	! should not be so small that it is impossible to fit all our photosynthesis
153	! levels in.
154
155	DO il=1,nlevels-1
156
157	ideal_level_thickness=(z_level(il+1)-z_level(il))&
158	/REAL(nl_photo,r_std)
159
160	DO ivm=1,nvm
161
162	IF(ideal_level_thickness .LE. min_level_sep(ivm))THEN
163
164	WRITE(numout,*) 'ERROR: We do not have enough space between the energy '
165	WRITE(numout,*) 'ERROR: levels to fit all the photosynthesis levels you '
166	WRITE(numout,*) 'ERROR: are asking for. Either reduce nl_photo or min_level_sep, '
167	WRITE(numout,*) 'ERROR: or increase the width of the levels for the energy budget.'
168
169	WRITE(numout,*) 'il, ivm: ',il,ivm
170	WRITE(numout,*) 'ideal_level_thickness, min_level_sep(ivm): ',&
171	ideal_level_thickness, min_level_sep(ivm)
172	WRITE(numout,*) 'z_level(il+1),z_level(il): ',z_level(il+1),z_level(il)
173	WRITE(numout,*) 'nl_photo: ',nl_photo
174
175	CALL ipslerr_p (3,'calculate_z_level_photo', 'Not enough space.','','')
176	ENDIF
177
178	ENDDO
179
180	ENDDO
181
182	DO ipts=1,npts
183	DO ivm=1,nvm
184
185	! Skip all of this if we have no leaf biomass
186	IF(biomass(ipts,ivm,ileaf,icarbon) .LE. min_stomate)THEN
187	z_level_photo(ipts,ivm,:)=zero
188	CYCLE
189	ENDIF
190
191	IF(is_tree(ivm))THEN
192
193	! We need to know the maximum height of the vegetation so that
194	! we don't calculate a lot of levels with nothing in them
195	! Notice that we have to use wood_to_height_eff here. The reason is that
196	! after coppicing, we may have zero aboveground wood, so our height calculated
197	! with wood_to_height would be zero.
198	heights(:)=wood_to_height_eff(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
199	max_height=MAXVAL(heights(:))
200
201	! Find the crown diameter and subtract this from the height,
202	! so that we know where the bottom of the canopy is.
203	cn_dia(:)=SQRT(4._r_std/pi*wood_to_cn_eff(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm))
204	cn_bottom(:)=heights(:)-cn_dia(:)
205	min_height=MINVAL(cn_bottom(:))
206
207	!IF(ld_laieff)THEN
208	IF((ivm == test_pft) .AND. (ipts == test_grid))THEN
209	WRITE(numout,*) 'This is the largest height of trees in z_level_photo'
210	WRITE(numout,*) 'max_height, min_height: ',max_height,min_height
211	WRITE(numout,*) 'heights: ',heights(:)
212	WRITE(numout,*) 'cn_dia: ',cn_dia(:)
213	! DO icir=1,ncirc
214	! WRITE(numout,*) 'cc_biomass: ',circ_class_biomass(ipts,ivm,icir,:,icarbon)
215	! ENDDO
216	ENDIF
217	!ENDIF
218
219	ELSE
220
221	! I took this from stomate_growth_fun_all
222	max_height = biomass_to_lai(biomass(ipts,ivm,ileaf,icarbon),ivm)&
223	* lai_to_height(ivm)
224
225	min_height = 0.0001_r_std ! arbitrary. Don't want it equal to zero because
226	! it won't enter the correct loop below.
227
228	IF(ld_laieff)THEN
229	IF((ivm == test_pft) .AND. (ipts == test_grid))THEN
230	WRITE(numout,*) 'This is the largest height of grasses/crops in z_level_photo'
231	WRITE(numout,*) 'max_height, min_height: ',max_height,min_height
232	WRITE(numout,*) 'lai_to_height: ',lai_to_height(ivm)
233	ENDIF
234	ENDIF
235
236	ENDIF
237
238	! Loop over the levels used by the energy budget. For each level, there are
239	! six possibilities, based on the location of the canopy vis-a-vis the location
240	! of this level height and the level above (if there is no level above, we take
241	! a boundary that is equal to a very big number).
242	! The word "slice" refers to the space between this level and the level above.
243	! 1) Only the bottom part of the canopy is in this slice.
244	! 2) Only the top part of the canopy is in this slice.
245	! 3) The whole canopy is below this slice.
246	! 4) The whole canopy is above this slice.
247	! 5) The whole canopy is inside this slice.
248	! 6) The top of the canopy is above the slice, the bottom is below the
249	! slice, and we just have the middle part in the slice.
250	! We adjust the photosynthesis levels in this slice to cover as much as the canopy as possible, assuming
251	! a mininum distince between the levels. This is done to avoid having a lot of levels with
252	! very little LAI, which could cause a problem in the photosyntheis routines.
253	! If there is too much LAI, the absorbed light saturates and we lose
254	! sensitivity in the photosynthesis routines. 0.5 seems to be a good LAI
255	! to have in a level, but this is just a feeling.
256
257	DO il=1,nlevels
258
259	lower_level_boundary=z_level(il)
260
261	! The lowest photosynthesis level should always be equal to
262	! the lower_level_boundary. In this way we make sure that
263	! we account for all the vegetation in this level. We need
264	! to make sure we don't overwrite this below. This means that
265	! the highest PS level in this level should always be at least
266	! min_level_sep from the upper_boundary, and that instead of
267	! nl_photo levels to play with, we only have nl_photo-1. This
268	! should not cause any divide by zero errors since there is
269	! a loop at the beginning of the subroutine to handle the trivial
270	! case where nl_photo=1.
271	z_level_photo(ipts,ivm,nl_photo*(il-1)+1)=lower_level_boundary
272
273	IF(il .NE. nlevels)THEN
274	upper_level_boundary=z_level(il+1)
275	ELSE
276	! This number is arbitary, but it just needs to be larger than
277	! the largest possible vegetation height. 10000 meters is a very
278	! tall tree.
279	upper_level_boundary=10000.0_r_std
280	ENDIF
281
282	IF((min_height .GT. lower_level_boundary .AND. &
283	min_height .LT. upper_level_boundary ) .AND. &
284	(max_height .GT. upper_level_boundary))THEN
285	! 1) This slice only has the bottom part of the canopy in it.
286	! Ideally, we would want the second lowest PS level to be just
287	! above the bottom of the canopy, and the highest PS level
288	! to be just below the upper boundary. If there is not
289	! enough of the canopy to do this and maintain our
290	! minimum level spacing, we just space the levels
291	! with this minimum starting from just below the upper boundary.
292
293	ideal_level_thickness = ((upper_level_boundary-min_level_sep(ivm))-&
294	(min_height+min_level_sep(ivm)))&
295	/REAL(nl_photo-2,r_std)
296
297	IF(ideal_level_thickness .GT. min_level_sep(ivm))THEN
298	! It worked, the levels are thick enough.
299	DO il_photo=nl_photo,2,-1
300	z_level_photo(ipts,ivm,nl_photo*(il-1)+il_photo)=&
301	(upper_level_boundary-min_level_sep(ivm))-&
302	ideal_level_thickness*REAL(nl_photo-il_photo,r_std)
303	ENDDO
304
305	ELSE
306
307	! Use the minimum level thickness starting from the
308	! upper boundary.
309	DO il_photo=nl_photo,2,-1
310	z_level_photo(ipts,ivm,nl_photo*(il-1)+il_photo)=&
311	(upper_level_boundary-min_level_sep(ivm))-&
312	min_level_sep(ivm)*REAL(nl_photo-il_photo,r_std)
313	ENDDO
314
315	ENDIF
316
317	ELSEIF(min_height .LT. lower_level_boundary .AND. &
318	((max_height .GT. lower_level_boundary) .AND. &
319	(max_height .LT. upper_level_boundary)))THEN
320	! 2) Only the top half of the canopy is in this slice. We need
321	! to do the same thing as case (1), but in the other direction.
322
323	ideal_level_thickness = ((max_height-min_level_sep(ivm))-&
324	(lower_level_boundary+min_level_sep(ivm)))&
325	/REAL(nl_photo-2,r_std)
326
327	IF(ideal_level_thickness .GT. min_level_sep(ivm))THEN
328	! It worked, the levels are thick enough.
329	DO il_photo=2,nl_photo
330	z_level_photo(ipts,ivm,nl_photo*(il-1)+il_photo)=&
331	(lower_level_boundary+min_level_sep(ivm))+&
332	ideal_level_thickness*REAL(il_photo-2,r_std)
333	ENDDO
334
335	ELSE
336
337	! Use the minimum level thickness starting from the
338	! upper boundary.
339	DO il_photo=2,nl_photo
340	z_level_photo(ipts,ivm,nl_photo*(il-1)+il_photo)=&
341	(lower_level_boundary+min_level_sep(ivm))+&
342	min_level_sep(ivm)*REAL(il_photo-2,r_std)
343	ENDDO
344
345	ENDIF
346
347	ELSEIF((min_height .GT. upper_level_boundary) .OR. &
348	(max_height .LT. lower_level_boundary) )THEN
349	! (3) and (4) The whole canopy is below the slice, or the whole canopy
350	! is above this slice. Whatever we do does not
351	! matter, since there is no LAI in the slice and therefore
352	! there will be no PS.
353
354	DO il_photo=2,nl_photo
355	z_level_photo(ipts,ivm,nl_photo*(il-1)+il_photo)=lower_level_boundary+&
356	min_level_sep(ivm)*REAL(il_photo-1,r_std)
357	ENDDO
358
359	ELSEIF((min_height .GT. lower_level_boundary) .AND. &
360	(max_height .LT. upper_level_boundary) )THEN
361	! (5) The whole canopy is inside this slice. This is the trickiest
362	! case, and unfortunately the most common. In fact, it's the only
363	! case which can exist with a single level energy budget.
364
365	ideal_level_thickness = ((max_height-min_level_sep(ivm))-&
366	(min_height+min_level_sep(ivm)))&
367	/REAL(nl_photo-2,r_std)
368
369	IF(ideal_level_thickness .GT. min_level_sep(ivm))THEN
370	! It worked, the levels are thick enough.
371	DO il_photo=2,nl_photo
372	z_level_photo(ipts,ivm,nl_photo*(il-1)+il_photo)=&
373	(min_height+min_level_sep(ivm))+&
374	ideal_level_thickness*REAL(il_photo-2,r_std)
375	ENDDO
376
377	ELSE
378
379	! Here we have a block of min_level_sep levels. Let's place the
380	! center of this block in the center of our canopy, and check to see
381	! if the ends poke past the boundaries of this slice.
382
383	block_thickness = REAL(nl_photo-2,r_std) * min_level_sep(ivm)
384	canopy_thickness = (max_height - min_height )
385	canopy_center = min_height + canopy_thickness / deux
386
387	IF( ((canopy_center + block_thickness/deux) .LT. &
388	(upper_level_boundary - min_level_sep(ivm))) .AND. &
389	(canopy_center - block_thickness/deux) .GT. &
390	(lower_level_boundary + min_level_sep(ivm)))THEN
391	! The block center can be placed at the center of the canopy, which means that
392	! the bottom of the block starts at canopy_center - block_thickness/2
393
394	DO il_photo=2,nl_photo
395	z_level_photo(ipts,ivm,nl_photo*(il-1)+il_photo)=&
396	canopy_center - block_thickness/deux+&
397	min_level_sep(ivm)*REAL(il_photo-2,r_std)
398	ENDDO
399
400	ELSEIF( (canopy_center - block_thickness/deux) .LE. &
401	(lower_level_boundary + min_level_sep(ivm)))THEN
402	! we have enough space up, but not down, so move the whole block upwards.
403	! All we have to do is start the block above the lower boundary level.
404
405	DO il_photo=2,nl_photo
406	z_level_photo(ipts,ivm,nl_photo*(il-1)+il_photo)=&
407	lower_level_boundary + &
408	min_level_sep(ivm)*REAL(il_photo-1,r_std)
409	ENDDO
410
411	ELSEIF( (canopy_center + block_thickness/deux) .GE. &
412	(upper_level_boundary - min_level_sep(ivm)))THEN
413	! we have enough space down, but not up, so move the whole block downwards
414	! All we need to do is start the block just below the upper boundary.
415
416	DO il_photo=nl_photo,2,-1
417	z_level_photo(ipts,ivm,nl_photo*(il-1)+il_photo)=&
418	upper_level_boundary - &
419	min_level_sep(ivm)*REAL(nl_photo-il_photo+1,r_std)
420	ENDDO
421
422	ELSE
423	! We should not be here!
424	WRITE(numout,*) 'ERROR: In calculate_z_photo_level, we have reached a '
425	WRITE(numout,*) 'ERROR: case that I have not foreseen!'
426	WRITE(numout,*) 'ipts,ivm,il: ',ipts,ivm,il
427	WRITE(numout,*) 'canopy_center: ',canopy_center
428	WRITE(numout,*) 'block_thickness: ',block_thickness
429	WRITE(numout,*) 'upper_level_boundary: ',upper_level_boundary
430	WRITE(numout,*) 'lower_level_boundary: ',lower_level_boundary
431	WRITE(numout,*) 'nl_photo: ',nl_photo
432
433	CALL ipslerr_p (3,'calculate_z_level_photo', 'Unforeseen case.','','')
434
435	ENDIF
436	ENDIF
437
438	ELSEIF((min_height .LT. lower_level_boundary) .AND. &
439	(max_height .GT. upper_level_boundary) )THEN
440	! (6) This slice is completely filled by the canopy. This is easy,
441	! we just evenly space the levels in this slice.
442
443	ideal_level_thickness = ((upper_level_boundary-min_level_sep(ivm))-&
444	(lower_level_boundary+min_level_sep(ivm)))&
445	/REAL(nl_photo-2,r_std)
446
447	DO il_photo=2,nl_photo
448	z_level_photo(ipts,ivm,nl_photo*(il-1)+il_photo)=&
449	(lower_level_boundary+min_level_sep(ivm))+&
450	ideal_level_thickness*REAL(il_photo-2,r_std)
451	ENDDO
452
453	ELSE
454	! Should not be here! This is a case I haven't thought of.
455	WRITE(numout,*) 'Should not be here in find photosynthesis levels!'
456	WRITE(numout,*) 'ipts,ivm,il: ',ipts,ivm,il
457	WRITE(numout,*) 'max_height: ',max_height
458	WRITE(numout,*) 'min_height: ',min_height
459	WRITE(numout,*) 'upper_level_boundary: ',upper_level_boundary
460	WRITE(numout,*) 'lower_level_boundary: ',lower_level_boundary
461	CALL ipslerr_p (3,'calculate_z_level_photo', 'Unforeseen case in ps levels.','','')
462
463	ENDIF ! All the possible locations of the canopy with respect to
464	! the level boundaries.
465
466	ENDDO
467
468	IF(ld_laieff .AND. ivm == test_pft .AND. ipts == test_grid)THEN
469	WRITE(numout,*) 'Printing the level heights used by photosynthesis.'
470	WRITE(numout,*) 'ipts,ivm,max_height: ',ipts,ivm,max_height
471	DO il=1,nlevels_tot
472	WRITE(numout,*) 'il,z_level_photo(ipts,ivm,il): ',&
473	il,z_level_photo(ipts,ivm,il)
474	ENDDO
475	ENDIF
476	ENDDO
477	ENDDO
478
479	IF (bavard.GE.2) WRITE(numout,*) 'Leaving calculate_z_level_photo'
480
481	END SUBROUTINE calculate_z_level_photo
482
483	!! ================================================================================================================================
484	!! SUBROUTINE : combine_lai_levels
485	!!
486	!>\BRIEF Collapses the LAI per photosynthesis level into an array
487	!! which only has the LAI per energy budget level.
488	!!
489	!! DESCRIPTION :
490	!!
491	!! RECENT CHANGE(S) : None
492	!!
493	!! MAIN OUTPUT VARIABLE(S): ::lai_per_level_temp
494	!!
495	!! REFERENCE(S) :
496	!!
497	!! FLOWCHART : None
498	!! \n
499	!_ ================================================================================================================================
500
501	SUBROUTINE combine_lai_levels(lai_per_level, lai_per_level_temp)
502
503	!! 0 Variable and parameter declaration
504
505	!! 0.1 Input variables
506
507	REAL(r_std),DIMENSION(:,:,:),INTENT(IN) :: lai_per_level !! The LAI per photosynthesis level
508	!! @tex $(m^{2} m^{-2})$
509
510
511	!! 0.2 Output variables
512	REAL(r_std),DIMENSION(:,:,:),INTENT(OUT) :: lai_per_level_temp !! The LAI per energy budget level
513	!! @tex $(m^{2} m^{-2})$
514
515	!! 0.3 Modified variables
516
517	!! 0.4 Local variables
518	INTEGER :: il, ilevel !! indices
519	!_ ================================================================================================================================
520
521	IF (bavard.GE.2) WRITE(numout,*) 'Entering combine_lai_levels'
522
523	lai_per_level_temp(:,:,:)=zero
524
525	DO ilevel=1,nlevels_tot
526	il=CEILING(REAL(ilevel)/REAL(nlevels_photo))
527	lai_per_level_temp(:,:,il)=lai_per_level_temp(:,:,il)+lai_per_level(:,:,ilevel)
528	ENDDO
529
530	IF (bavard.GE.2) WRITE(numout,*) 'Leaving combine_lai_levels'
531
532	END SUBROUTINE combine_lai_levels
533
534	!! ================================================================================================================================
535	!! SUBROUTINE : find_lai_per_level
536	!!
537	!>\BRIEF Calculates the LAI per vertical level of the canopy.
538	!!
539	!! DESCRIPTION : Given a total biomass, number of individuals, size of individuals,
540	!! and set of vertical layer boundaries, this routine calculates the
541	!! the leaf area index (LAI) found between the layer boundaries. It
542	!! first makes a call to the stand structure in order to find out some
543	!! canopy properties like height and canopy diameter, and then uses geometry
544	!! to determine the amount of LAI present.
545	!!
546	!! RECENT CHANGE(S) : None
547	!!
548	!! MAIN OUTPUT VARIABLE(S): ::lai_per_level
549	!!
550	!! REFERENCE(S) :
551	!!
552	!! FLOWCHART : None
553	!! \n
554	!_ ================================================================================================================================
555
556	SUBROUTINE find_lai_per_level(npts, nlevels_loc, z_level_photo, &
557	biomass, circ_class_biomass, circ_class_n, lai_per_level, max_height_store)
558
559	!! 0 Variable and parameter declaration
560
561	!! 0.1 Input variables
562
563	INTEGER(i_std),INTENT(IN) :: npts !! Domain size - number of pixels (unitless)
564	INTEGER(i_std),INTENT(IN) :: nlevels_loc !! The number of levels we divide the LAI over (unitless)
565	REAL(r_std),DIMENSION(:,:,:),INTENT(IN) :: z_level_photo !! The height of the bottom of each of our levels @tex $(m)$ @endtex
566	!! Note that level 1 is closest to the ground
567	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(IN) :: circ_class_biomass !! Biomass of the componets of the model
568	!! tree within a circumference
569	!! class @tex $(gC ind^{-1})$ @endtex
570	REAL(r_std), DIMENSION(:,:,:), INTENT(IN) :: circ_class_n !! Number of trees within each circ class
571	!! class @tex $(m^{-2})$ @endtex
572	REAL(r_std), DIMENSION(:,:,:,:), INTENT(IN) :: biomass !! Stand level biomass @tex $(gC m^{-2})$ @endtex
573
574	!! 0.2 Output variables
575	REAL(r_std), DIMENSION(:,:,:), INTENT(OUT) :: lai_per_level !! This is the LAI per vertical level
576	!! @tex $(m^{2} m^{-2})$ @endtex
577	REAL(r_std), DIMENSION(npts, nvm), INTENT(OUT) :: max_height_store
578
579
580	!! 0.3 Modified variables
581
582	!! 0.4 Local variables
583	INTEGER :: ipts, ivm, il, il_photo, &
584	nlevels_loc_photo, icir !! indices
585	REAL(r_std) :: max_height !! The height of the tallest tree
586	!! @tex $(m)$ @endtex
587	REAL(r_std) :: ind_loc !! The number of individuals
588	!! @tex $(m^{-2})$ @endtex
589	REAL(r_std) :: upper_boundary !! The top of the current level
590	!! @tex $(m)$ @endtex
591	REAL(r_std) :: lower_boundary !! The bottom of the current level
592	!! @tex $(m)$ @endtex
593	REAL(r_std) :: lai_temp !! LAI of the current PFT/grid square
594	!! @tex $(m^{2} m^{-2})$ @endtex
595	REAL(r_std),DIMENSION(npts,nvm,ncirc,ndist_types) :: values !! An array which holds various canopy parameters
596	REAL(r_std),DIMENSION(npts,nvm,nlevels_loc) :: PartialCrownVolumeL !! The fraction of the crown volume found
597	!! in this level averaged over all circ classes
598	!! @tex $(m^{3})$ @endtex
599	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) :: PartialCrownVolume_h !! The fraction of the crown volume found
600	!! in this level for each circ class
601	!! @tex $(m^{3})$ @endtex
602	REAL(r_std),DIMENSION(npts,nvm) :: CrownVolumeL !! The average crown volume of the stand.
603	!! @tex $(m^3)$ @endtex
604	REAL(r_std),DIMENSION(npts,nvm) :: FAVD !! Foliage area volume density
605	!! @tex $(m^2 (leaf) m^{-3})$ @endtex
606
607	!_ ================================================================================================================================
608
609	IF (bavard.GE.2) WRITE(numout,*) 'Entering find_lai_per_level'
610
611	CALL derive_biomass_quantities(npts, nvm, ncirc, circ_class_n, &
612	circ_class_biomass, values)
613
614	PartialCrownVolumeL(:,:,:)=zero
615	PartialCrownVolume_h(:,:,:,:)=zero
616	CrownVolumeL(:,:)=zero
617
618	DO ipts=1,npts
619	DO ivm=1,nvm
620
621	! Skip all of this if we have no leaf biomass
622	IF(biomass(ipts,ivm,ileaf,icarbon) .LE. min_stomate)THEN
623	lai_per_level(ipts,ivm,:)=zero
624	CYCLE
625	ENDIF
626
627	ind_loc=SUM(circ_class_n(ipts,ivm,:))
628
629	IF(ind_loc .LT. min_stomate)THEN
630	WRITE(numout,*) 'How do we have biomass but no individuals?'
631	WRITE(numout,*) 'ipts,ivm: ',ipts,ivm
632	WRITE(numout,*) 'biomass(ipts,ivm,ileaf,icarbon),ind_loc: ',&
633	biomass(ipts,ivm,ileaf,icarbon),ind_loc
634	CALL ipslerr_p (3,'find_lai_per_level', 'Biomass but no individuals.','','')
635	ENDIF
636
637	CALL calculate_favd(ipts,ivm, values(ipts,ivm,:,icnvol), circ_class_biomass, &
638	circ_class_n, biomass, favd)
639
640	! Find the LAI of this PFT.
641	lai_temp = biomass_to_lai(biomass(ipts,ivm,ileaf,icarbon), ivm)
642
643
644	! This is done differently for forests and grass/crops, since grass
645	! and crops will not have a canopy volume. We just assume a solid
646	! block of vegetation with a given height.
647
648	IF(is_tree(ivm))THEN
649
650	IF(ld_laieff .AND. test_pft == ivm .AND. ipts == test_grid)&
651	WRITE(6,*) 'On tree ',ivm,lai_temp,favd(ipts,ivm)
652
653	max_height = MAXVAL(values(ipts,ivm,:,iheight))
654
655	! We need to integrate over the height distribution. This is just a
656	! simple rectangular integration, and it gives us the average crown
657	! volume for the whole stand.
658	DO icir=1,ncirc
659	CrownVolumeL(ipts,ivm)=CrownVolumeL(ipts,ivm)+&
660	values(ipts,ivm,icir,icnvol)*&
661	circ_class_n(ipts,ivm,icir)/ind_loc
662	IF(test_pft == ivm .AND. ld_laieff .AND. ipts == test_grid)THEN
663	WRITE(numout,*) 'Volume for icir: ',icir,values(ipts,ivm,icir,icnvol)
664	ENDIF
665	ENDDO
666
667	! Now I need to find the average crown volume per layer. This is not the
668	! same thing as taking the average crown volume and dividing it up
669	! over the layers, since the crowns are at different heights.
670
671	DO il=1,nlevels_loc
672
673	IF(il .NE. nlevels_loc)THEN
674	upper_boundary=z_level_photo(ipts,ivm,il+1)
675	ELSE
676	! This height is arbitrary. It just has to include the whole canopy.
677	upper_boundary=z_level_photo(ipts,ivm,nlevels_loc)+max_height
678	ENDIF
679	lower_boundary=z_level_photo(ipts,ivm,il)
680
681	! For this level, how much of each canopy is here? We will also
682	! integrate this value over the whole height distribution, thus
683	! finding the average canopy volume in this layer.
684	DO icir=1,ncirc
685
686	PartialCrownVolume_h(ipts,ivm,icir,il)=&
687	partial_spheroid_vol(upper_boundary, lower_boundary, &
688	values(ipts,ivm,icir,icndiaver)/deux, &
689	values(ipts,ivm,icir,icndiahor)/deux, &
690	values(ipts,ivm,icir,iheight))
691	IF(ld_laieff .AND. test_pft == ivm .AND. test_grid == ipts)THEN
692	WRITE(numout,*) 'ipts,ivm,icir,il: ',ipts,ivm,icir,il
693	WRITE(numout,*) 'PartialCrownVolume_h(ipts,ivm,icir,il): ',&
694	PartialCrownVolume_h(ipts,ivm,icir,il)
695	WRITE(numout,*) 'crown diameter: ', values(ipts,ivm,icir,icndiaver)
696	WRITE(numout,*) 'tree height: ', values(ipts,ivm,icir,iheight)
697	WRITE(numout,*) 'circ_class_n(ipts,ivm,icir)',circ_class_n(ipts,ivm,icir)
698	ENDIF
699
700	! It's important to note here that the units of PartialCrownVolume_h
701	! are m**3/tree, since the canopy properties that we passed to
702	! partial_spheroid_vol are also per tree. We need to get it in
703	! m3/m2 (land area).
704
705	PartialCrownVolumeL(ipts,ivm,il)=PartialCrownVolumeL(ipts,ivm,il)+&
706	PartialCrownVolume_h(ipts,ivm,icir,il)*&
707	circ_class_n(ipts,ivm,icir)
708
709	END DO
710
711	! Now I have the average crown volume per level. I also have
712	! the foliage area volume density. The amount of LAI in this
713	! level is just the combination of the two.
714	lai_per_level(ipts,ivm,il)=PartialCrownVolumeL(ipts,ivm,il)*&
715	FAVD(ipts,ivm)
716
717	ENDDO
718
719	IF(ld_laieff .AND. ipts == test_grid .AND. ivm == test_pft)THEN
720	WRITE(numout,*) 'CrownVolumeL: ',CrownVolumeL(ipts,ivm)
721	ENDIF
722
723	ELSE
724
725	! For grasses and crops.
726
727	! This is not ideal, since we already had to calculate lai_temp
728	! and max_height to find the FAVD. If there is a speed problem,
729	! we can change it.
730
731	! What is the height of our grass/crop?
732
733	! This converts LAI to a height. I took this from functional_allocation.
734	max_height = lai_temp * lai_to_height(ivm)
735
736
737	! Find out how much of this vegetation cube is in this layer
738	! These calculations here include some implict factors of
739	! m2/m2, so be aware of that when calculating the levels.
740	! All heights are multiplied by 1x1 m to find a volume, but
741	! that is not explictly shown.
742	DO il=1,nlevels_loc
743
744	IF(il .NE. nlevels_loc)THEN
745	upper_boundary=z_level_photo(ipts,ivm,il+1)
746	ELSE
747	! This height is arbitrary. It just has to include the whole canopy.
748	upper_boundary=z_level_photo(ipts,ivm,nlevels_loc)+max_height
749	ENDIF
750	lower_boundary=z_level_photo(ipts,ivm,il)
751
752	IF(max_height .GE. upper_boundary)THEN
753	! This whole level is full of vegetation.
754	lai_per_level(ipts,ivm,il)=&
755	(upper_boundary-lower_boundary)*FAVD(ipts,ivm)
756
757	ELSEIF(max_height .LE. lower_boundary)THEN
758
759	! There is no vegetation in this level.
760	lai_per_level(ipts,ivm,il)=zero
761
762	ELSE
763
764	! This level is partially full.
765	lai_per_level(ipts,ivm,il)=&
766	(max_height-lower_boundary)*FAVD(ipts,ivm)
767
768	ENDIF ! checking where max_height is
769
770	ENDDO
771
772	ENDIF
773
774	IF(test_pft == ivm .AND. test_grid == ipts .AND. ld_laieff)THEN
775	WRITE(numout,*) 'LAI per level: ',ipts,ivm
776	DO il=nlevels_loc,1,-1
777	WRITE(numout,*) 'il,lai_per_level(ipts,ivm,il): ',&
778	il,lai_per_level(ipts,ivm,il)
779	ENDDO
780	ENDIF
781
782	! Test to see if our LAI is equal to the sum over all the levels.
783	IF( ABS(lai_temp - SUM(lai_per_level(ipts,ivm,:))) &
784	.GT. 1e-5)THEN
785	WRITE(numout,) '!****************************************************'
786	WRITE(numout,*) 'Problem in LAI levels! '
787	WRITE(numout,*) 'ipts,ivm: ',ipts,ivm
788	WRITE(numout,*) 'lai_temp, SUM(lai_per_level(ipts,ivm,:)): ',&
789	lai_temp,SUM(lai_per_level(ipts,ivm,:))
790	DO il=nlevels_loc,1,-1
791	WRITE(numout,*) 'il,lai_per_level(ipts,ivm,il): ',&
792	il,lai_per_level(ipts,ivm,il)
793	ENDDO
794	! There can be issues using the functions in the print statements. If a print
795	! statement in a function is activated when the function call itself is in a print
796	! statement, there is a crash.
797	WRITE(numout,*) 'Check whether canopy dimensions exceed the height '
798	! WRITE(numout,*) 'diameter, ', wood_to_dia(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
799	WRITE(numout,*) 'diameter, ', values(ipts,ivm,:,idiameter)
800	! WRITE(numout,*) 'height: ', wood_to_height(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
801	WRITE(numout,*) 'height: ', values(ipts,ivm,:,iheight)
802	! WRITE(numout,*) 'crown area: ', wood_to_cn(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
803	WRITE(numout,*) 'crown area: ', values(ipts,ivm,:,icnarea)
804	! WRITE(numout,) 'crown diameter: ', (4./piwood_to_cn(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm))**(1./2.)
805	WRITE(numout,*) 'crown diameter - ver: ', values(ipts,ivm,:,icndiaver)
806	WRITE(numout,*) 'crown diameter - hor: ', values(ipts,ivm,:,icndiahor)
807	WRITE(numout,*) 'crown volume: ', values(ipts,ivm,:,icnvol)
808	WRITE(numout,) 'crown diameter used in this routine: ', (6./pivalues(ipts,ivm,:,icnvol))**(1./3.)
809	WRITE(numout,) '!****************************************************'
810	CALL ipslerr_p (2,'find_lai_per_level', 'Inconsistency with LAI in levels.','','')
811	ENDIF ! check to see if the LAI per level is consistent with the total LAI
812
813	max_height_store(ipts, ivm) = max_height
814
815
816	ENDDO ! loop over PFTs
817	ENDDO ! loop over points
818
819	IF (bavard.GE.2) WRITE(numout,*) 'Leaving find_lai_per_level'
820
821	END SUBROUTINE find_lai_per_level
822
823
824	!! ================================================================================================================================
825	!! SUBROUTINE : fitting_laieff
826	!!
827	!>\BRIEF Fits the effective LAI as a function of the solar angle
828	!!
829	!! DESCRIPTION : Calculation of the effective LAI is a very expensive part of the
830	!! model, due to the loops over points, PFTs, circ classes, and levels,
831	!! in addition to using trig functions. We need this value at
832	!! every time step (up to 48 times a day). Thankfully, the function
833	!! doesn't change a lot, in particular when the sun is high in the sky.
834	!! Therefore, we will calculate a few points exactly, and then fit
835	!! a function to these points to predict the value of the LAIeff at any
836	!! solar angle. This method seems to work well for low solar zenith
837	!! angles, which are the most important since they will have the most
838	!! incoming solar radiation.
839	!!
840	!! NOTE:
841	!!
842	!! RECENT CHANGE(S) : None
843	!!
844	!! MAIN OUTPUT VARIABLE(S): ::laieff_fit
845	!!
846	!! REFERENCE(S) :
847	!!
848	!! FLOWCHART : None
849	!! \n
850	!_ ================================================================================================================================
851
852	SUBROUTINE fitting_laieff(npts, nlevels_loc, z_level_photo, circ_class_biomass, &
853	circ_class_n, veget_max, biomass, lai_per_level, laieff_fit, h_array_out, z_array)
854
855	!! 0 Variable and parameter declaration
856
857	!! 0.1 Input variables
858
859	INTEGER(i_std),INTENT(IN) :: npts !! Domain size - number of pixels (unitless)
860	INTEGER(i_std),INTENT(IN) :: nlevels_loc !! The number of vertical levels
861	REAL(r_std),DIMENSION(:,:,:),INTENT(IN) :: z_level_photo !! the height of the bottom of each of our levels
862	!! @tex $(m)$ @endtex
863	!! note that level 1 is closest to the ground
864	REAL(r_std), DIMENSION(:,:,:,:,:), &
865	INTENT(IN) :: circ_class_biomass !! Biomass of the componets of the model
866	!! tree within a circumference class
867	!! @tex $(gC ind^{-1})$ @endtex
868	REAL(r_std), DIMENSION(:,:,:), INTENT(IN) :: circ_class_n !! Number of trees within each circ class
869	!! @tex $(m^{-2})$ @endtex
870	REAL(r_std),DIMENSION(:,:),INTENT(IN) :: veget_max !! Maximum fraction of vegetation type including
871	!! non-biological fraction (unitless)
872	REAL(r_std), DIMENSION(:,:,:,:), INTENT(IN) :: biomass !! Stand level biomass
873	!! @tex $(gC m^{-2})$ @endtex
874	REAL(r_std), DIMENSION(:,:,:), INTENT(IN) :: lai_per_level !! This is the LAI per vertical level
875	!! @tex $(m^{2} m^{-2})$
876	!! 0.2 Output variables
877	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(OUT) :: laieff_fit !! Fitted parameters for the effective LAI
878
879
880	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc), INTENT(OUT) &
881	:: h_array_out !! An output of h_array, to use in sechiba
882	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc), INTENT(OUT) &
883	:: z_array !! Height above soil of the Pgap points.
884	!! @tex $(m)$ @endtex
885
886
887	!! 0.3 Modified variables
888
889	!! 0.4 Local variables
890
891	! Idealy, the next two values would be externalised.
892	INTEGER,PARAMETER :: nangles=8 ! the number of angles to fit the parameterization of
893	! the effective LAI.
894	! REAL(r_std),DIMENSION(nangles),PARAMETER :: param_angles = (/ 0.1, 30.0, 60.0, 89.9 /) ! the angles to
895	REAL(r_std),DIMENSION(nangles) :: param_angles ! the angles to
896	! REAL(r_std),DIMENSION(nangles) :: param_angles
897	! parameterize for...EXTERNALIZE
898
899	REAL(r_std),DIMENSION(nlevels_tot,npts,nvm,nangles) :: laieff !! Leaf Area Index Effective converts 3D lai into
900	!! 1D lai for two stream radiation transfer model
901	!! @tex $(m^{2} m^{-2})$ @endtex
902	REAL(r_std), DIMENSION(npts) :: solar_angle !! The solar zenith angle of every grid square
903	!! @tex $(rad)$ @endtex
904	INTEGER(i_std) :: ilevel, ipts, icir, ivm, &
905	iangle, jangle !! index (unitless)
906
907	INTEGER(i_std) :: nzero !! Counter to keep track of the number of
908	!! LAIeff with value zero
909	INTEGER(i_std) :: temp_index !! An index which allows us to fit only non-zero
910	!! values
911	REAL(r_std),DIMENSION(npts) :: test_angles !! Converts the param_angles (degrees) into radians
912	!! @tex $(rad)$ @endtex
913	REAL(r_std),DIMENSION(npts,nvm,nlevels_loc) &
914	:: jz_array2 !! Same as z_array, but one less dimension.
915	!! @tex $(m)$ @endtex (local because the call to effective_lai requires it)
916
917	!_ ================================================================================================================================
918
919	IF (bavard.GE.2) WRITE(numout,*) 'Entering fitting_laieff'
920
921	!! 1. Initialize check for surface area conservation
922	! Veget_max is a INTENT(in) variable and can therefore
923	! not be changed during the course of this subroutine
924	! No need to check whether the subroutine preserves the
925	! total surface area of the pixel.
926
927	!! 2. Calculations
928	! We want to parameterize the effective LAI as a function of
929	! the solar zenith angle for every level, grid square, and
930	! PFT type. Let's choose a series of angles to use.
931	! There are somtimes numerical issues with using exactly
932	! 90 degrees or exactly 0 degrees.
933
934	! param_angles = (/ 5.0, 15.0, 25.0, 35.0, 45.0, 55.0, 65.0, 75.0, 85.0 /)
935	! Using high SZA (near the horizon) seems to cause problems
936	! with the fitting, as strange behavior happens there. Since there
937	! is not a lot of light, let's use smaller values.
938	param_angles = (/ 5.0, 15.0, 25.0, 35.0, 45.0, 55.0, 65.0, 75.0 /)
939
940
941	! Right now, I'm doing this the brute force way. It will probably
942	! need to be refined later.
943
944	DO iangle=1,nangles
945	test_angles(:)=COS(param_angles(iangle)/180.0*pi)
946	IF(ld_laieff)&
947	WRITE(numout,*) 'Calculating for test angle: ',param_angles(iangle)
948	CALL effective_lai(npts, nlevels_loc, z_level_photo, circ_class_biomass, &
949	circ_class_n, veget_max, test_angles, biomass, lai_per_level, &
950	laieff(:,:,:,iangle), jz_array2, h_array_out, z_array)
951	IF(ld_laieff)THEN
952	DO ipts=1,npts
953	DO ivm=1,nvm
954	IF(ipts == test_grid .AND. ivm == test_pft)THEN
955	WRITE(numout,*) 'LAIEFF for angle ',param_angles(iangle)
956	DO ilevel=nlevels_tot,1,-1
957	WRITE(numout,*) ilevel,laieff(ilevel,ipts,ivm,iangle)
958	ENDDO
959	ENDIF
960	ENDDO
961	ENDDO
962	ENDIF
963	ENDDO
964
965	param_angles(:)=COS(param_angles(:)/180.0_r_std*pi)
966
967	! Now we have all the data points, so we can fit the function.
968	DO ipts=1,npts
969	DO ivm=1,nvm
970	DO ilevel=1,nlevels_loc
971
972	! The 0.001 is completely arbitrary here.
973	IF(SUM(laieff(ilevel,ipts,ivm,:)) .GT. 0.001_r_std )THEN
974
975	! We have a choice to make. Sometimes the laieff for highest zenith angles
976	! will be zero, even if the previous points were not. This seems to happen
977	! when no light gets through, which is normal for lower vegetation
978	! levels as the sun gets closer to the horizon. The pattern of points is
979	! something like 1,1,1,1,1,1.9,4.0,0.0. This causes a problem for
980	! fitting, and we don't care too much about this largest angle. If we
981	! can fit the first points well, and the last point is greater than 4.0,
982	! we will probably be okay since there is not a lot of sunlight wen the
983	! sun is near the horizon in any case. So let's run a test.
984	IF(laieff(ilevel,ipts,ivm,nangles) .LT. min_stomate)THEN
985
986	! We know that the sum of all elements is greater than our cutoff
987	! above, so we need to find the final substantial element.
988	! Count backwards to see how many zero elements we have
989	! at the end.
990	nzero=1
991	DO iangle=nangles-1,1,-1
992	IF(laieff(ilevel,ipts,ivm,iangle) .LT. min_stomate)THEN
993	nzero=nzero+1
994	ELSE
995	EXIT
996	ENDIF
997	ENDDO
998
999	IF(nzero .GE. nangles)THEN
1000	WRITE(numout,*) 'ERROR: Something went wrong in fitting the laieff!'
1001	WRITE(numout,*) 'ipts,ivm,ilevel: ',ipts,ivm,ilevel
1002	WRITE(numout,*) 'nzero,nangles: ',nzero,nangles
1003	WRITE(numout,*) 'laieff(ilevel,ipts,ivm,:): ',laieff(ilevel,ipts,ivm,:)
1004	CALL ipslerr_p (3,'fitting_laieff', 'All zero values in the fit.','','')
1005	ENDIF
1006
1007	! This three is also arbitrary.
1008	IF(nzero .GE. 3)THEN
1009	WRITE(numout,*) 'ERROR: A surprising number of zero values in laieff!'
1010	WRITE(numout,*) 'ipts,ivm,ilevel: ',ipts,ivm,ilevel
1011	WRITE(numout,*) 'nzero,nangles: ',nzero,nangles
1012	WRITE(numout,*) 'laieff(ilevel,ipts,ivm,:): ',laieff(ilevel,ipts,ivm,:)
1013	CALL ipslerr_p (3,'fitting_laieff', 'Too many zero values in the fit.','','')
1014	ENDIF
1015
1016	! Only use the first non-zero values for the fitting
1017	temp_index=nangles-nzero
1018	CALL fit_laieff_to_points(temp_index,param_angles(1:temp_index),&
1019	laieff(ilevel,ipts,ivm,1:temp_index), &
1020	laieff_fit(ipts,ivm,ilevel),ipts,ivm,ilevel)
1021	ELSE
1022	CALL fit_laieff_to_points(nangles,param_angles,laieff(ilevel,ipts,ivm,:), &
1023	laieff_fit(ipts,ivm,ilevel),ipts,ivm,ilevel)
1024	ENDIF
1025
1026	ELSE
1027
1028	! Not enough vegetation to fit it. Everything is zero.
1029	laieff_fit(ipts,ivm,ilevel)%a=zero
1030	laieff_fit(ipts,ivm,ilevel)%b=zero
1031	laieff_fit(ipts,ivm,ilevel)%c=zero
1032	laieff_fit(ipts,ivm,ilevel)%d=zero
1033	laieff_fit(ipts,ivm,ilevel)%e=zero
1034	ENDIF
1035	END DO
1036	ENDDO
1037	ENDDO
1038
1039	IF (bavard.GE.2) WRITE(numout,*) 'Leaving fitting_laieff'
1040
1041	END SUBROUTINE fitting_laieff
1042
1043
1044	!! ================================================================================================================================
1045	!! SUBROUTINE : fit_laieff_to_points
1046	!!
1047	!>\BRIEF
1048	!!
1049	!! DESCRIPTION :
1050	!!
1051	!! NOTE:
1052	!!
1053	!! RECENT CHANGE(S) : None
1054	!!
1055	!! MAIN OUTPUT VARIABLE(S): ::laieff_fit
1056	!!
1057	!! REFERENCE(S) :
1058	!!
1059	!! FLOWCHART : None
1060	!! \n
1061	!_ ================================================================================================================================
1062
1063	SUBROUTINE fit_laieff_to_points(nangles, test_angles, laieff, laieff_fit, ipts, ivm, ilevel)
1064
1065	!! 0 Variable and parameter declaration
1066
1067	!! 0.1 Input variables
1068
1069	INTEGER(i_std),INTENT(IN) :: nangles !! The number of angles (data points)
1070	INTEGER(i_std),INTENT(IN) :: ipts !! The grid square
1071	INTEGER(i_std),INTENT(IN) :: ivm !! The PFT number
1072	INTEGER(i_std),INTENT(IN) :: ilevel !! The canopy level
1073	REAL(r_std),DIMENSION(:),INTENT(IN) :: test_angles !! the cosine of the solar angles
1074	REAL(r_std),DIMENSION(:),INTENT(IN) :: laieff !! the laieff at each of these angles
1075
1076	!! 0.2 Output variables
1077	TYPE(laieff_type),INTENT(out) :: laieff_fit !! Fitted parameters for the effective LAI
1078
1079	!! 0.3 Modified variables
1080
1081	!! 0.4 Local variables
1082	REAL(r_std) :: x,y,one,x2,x3,x4,x5,x6,x7,x8,xy,x2y,x3y,x4y
1083	INTEGER :: iangle
1084	REAL(r_std),DIMENSION(5) :: vector_b
1085	REAL(r_std),DIMENSION(5,5) :: matrix_a
1086	REAL(r_std) :: mean,sse,sst,r_corr, max_diff, diff
1087
1088	!_ ================================================================================================================================
1089
1090	IF (bavard.GE.2) WRITE(numout,*) 'Entering fit_laieff_to_points'
1091
1092	! We have a set of angles, test_angles. We also have the
1093	! effective LAI for each of these angles. The structure
1094	! laieff_fit gives the parameters and functional form.
1095	! If you wish to change the functional form, you need to
1096	! change this routine, the lai_type, and the function
1097	! below which evaluates the LAIeff as a function of
1098	! the angle.
1099
1100	! We are going to use a least squares method to fit a parabola
1101	! of the form a+bx+cx*2+dx*3+ex**4. This means we will need to solve
1102	! a system of linear equations, A x = b.
1103	! In this notation, y is the sum over all angles for the laieff,
1104	! x is the sum for the angles, x2 is the sum of the square of the
1105	! angles, etc.
1106
1107	y=zero
1108	xy=zero
1109	x2y=zero
1110	x3y=zero
1111	x4y=zero
1112	one=zero
1113	x=zero
1114	x2=zero
1115	x3=zero
1116	x4=zero
1117	x5=zero
1118	x6=zero
1119	x7=zero
1120	x8=zero
1121	DO iangle=1,nangles
1122	y=y+laieff(iangle)
1123	xy=xy+test_angles(iangle)*laieff(iangle)
1124	x2y=x2y+test_angles(iangle)*2laieff(iangle)
1125	x3y=x3y+test_angles(iangle)*3laieff(iangle)
1126	x4y=x4y+test_angles(iangle)*4laieff(iangle)
1127	x=x+test_angles(iangle)
1128	x2=x2+test_angles(iangle)**2
1129	x3=x3+test_angles(iangle)**3
1130	x4=x4+test_angles(iangle)**4
1131	x5=x5+test_angles(iangle)**5
1132	x6=x6+test_angles(iangle)**6
1133	x7=x7+test_angles(iangle)**7
1134	x8=x8+test_angles(iangle)**8
1135	ENDDO
1136	one=REAL(nangles)
1137
1138	! Now we create our system of equations.
1139	vector_b(1:5) = (/ y, xy, x2y, x3y, x4y /)
1140	! From what I see in the gauss_jordan routine, the format is
1141	! matrix(row, column), which matches standard notation. In this
1142	! case the matrix is symmetric, anyway.
1143	matrix_a(1,1:5) = (/ one, x, x2, x3, x4 /)
1144	matrix_a(2,1:5) = (/ x, x2, x3, x4, x5 /)
1145	matrix_a(3,1:5) = (/ x2, x3, x4, x5, x6 /)
1146	matrix_a(4,1:5) = (/ x3, x4, x5, x6, x7 /)
1147	matrix_a(5,1:5) = (/ x4, x5, x6, x7, x8 /)
1148
1149	! TEST... should give -2, 3, 1...it does
1150	! matrix_a(1,1:3) = (/ -3, 2, -6 /)
1151	! matrix_a(2,1:3) = (/ 5, 7, -5 /)
1152	! matrix_a(3,1:3) = (/ 1, 4, -2 /)
1153	! vector_b(1:3) = (/ 6, 6, 8 /)
1154
1155	CALL gauss_jordan_method(SIZE(vector_b),matrix_a,vector_b)
1156
1157	! These are our parameters.
1158	laieff_fit%a=vector_b(1)
1159	laieff_fit%b=vector_b(2)
1160	laieff_fit%c=vector_b(3)
1161	laieff_fit%d=vector_b(4)
1162	laieff_fit%e=vector_b(5)
1163
1164	!++++++ CHECK ++++++++
1165	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1166	! These are some checks that we do to make sure our fitting
1167	! is decent. I'm going to leave these on during parameterization,
1168	! but we might want to remove it during the real runs. Things
1169	! should be tested well before then, and it will only be freak
1170	! chance that causes it to fail afterwards. If this happens
1171	! on one day in one hundred years, we don't care so much.
1172
1173	! Calculate a regression coefficienct
1174	sst=zero
1175	sse=zero
1176	mean=zero
1177	DO iangle=1,nangles
1178	mean=mean+laieff(iangle)
1179	ENDDO
1180	mean=mean/REAL(nangles)
1181	max_diff=zero
1182	DO iangle=1,nangles
1183	sst=sst+(laieff(iangle)-mean)**2
1184	diff=laieff(iangle)-calculate_laieff_fit(test_angles(iangle),laieff_fit)
1185	sse=sse+diff**2
1186	IF(diff .GT. max_diff) max_diff=diff
1187	ENDDO
1188
1189	IF(sse .GT. sst)THEN
1190	r_corr = zero
1191	ELSEIF(sst .EQ. zero)THEN
1192	r_corr = zero
1193	ELSE
1194	r_corr=SQRT(un-sse/sst)
1195	ENDIF
1196	IF(ld_laieff .AND. (ipts == test_grid) .AND. (ivm == test_pft)) &
1197	WRITE(numout,*) 'Regression coefficient: ',r_corr
1198
1199	IF(r_corr .LT. 0.95)THEN
1200
1201	! Sometimes the correlation coefficient is not great, but the difference
1202	! in the values is really small. We don't care so much if our values
1203	! are off by 0.001.
1204	IF(ld_laieff .AND. (ipts == test_grid) .AND. (ivm == test_pft)) &
1205	WRITE(numout,*) 'max_diff ',max_diff
1206	IF(max_diff .GT. 0.01)THEN
1207	WRITE(numout,*) 'Difference is too big for fitting!'
1208	WRITE(numout,*) 'max_diff ',max_diff
1209	WRITE(numout,*) 'r_corr ',r_corr
1210	WRITE(numout,*) 'ipts,ivm,ilevel: ',ipts,ivm,ilevel
1211	WRITE(numout,*) 'COS(angle) Real value Fitted value'
1212	DO iangle=1,nangles
1213	WRITE(1600,'(10F20.10)') test_angles(iangle),&
1214	laieff(iangle),calculate_laieff_fit(test_angles(iangle),laieff_fit)
1215	WRITE(numout,'(10F20.10)') test_angles(iangle),&
1216	laieff(iangle),calculate_laieff_fit(test_angles(iangle),laieff_fit)
1217	ENDDO
1218	CALL ipslerr_p (3,'fitting_laieff', 'Difference is too big for fitting.','','')
1219	ENDIF
1220	ENDIF
1221	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1222
1223	! Some debugging information
1224	IF(ld_laieff .AND. (ivm == test_pft) .AND. (ipts == test_grid))THEN
1225	! Let's see how closely they match
1226	DO iangle=1,nangles
1227	IF(SUM(laieff(:)) .GT. 0.1 )THEN
1228	WRITE(numout,*) 'TESTING fitting LAIEFF ',ipts,ivm,ilevel
1229	WRITE(numout,'(10F16.8)') test_angles(iangle),&
1230	laieff(iangle),calculate_laieff_fit(test_angles(iangle),laieff_fit)
1231	ENDIF
1232	ENDDO
1233	ENDIF
1234
1235	IF (bavard.GE.2) WRITE(numout,*) 'Leaving fit_laieff_to_points'
1236
1237	END SUBROUTINE fit_laieff_to_points
1238
1239	!! ================================================================================================================================
1240	!! SUBROUTINE : effective_lai
1241	!!
1242	!>\BRIEF Calculate effective lai for the 1-D two stream
1243	!! radiative canopy transfer model
1244	!!
1245	!! DESCRIPTION : The effective LAI is computed to be used with the two-stream
1246	!! albedo model described by Pinty 2006, by using an equation in Pinty 2004
1247	!! and the Pgap method of Haverd et al...code was taken from the latter group
1248	!! and cleaned up and reformatted...these routines were valided against the original
1249	!! code and gave identical results for Pgap (the original code was single precision, while
1250	!! these results are double precision):
1251	!!
1252	!! NOTE:
1253	!!
1254	!! RECENT CHANGE(S) : None
1255	!!
1256	!! MAIN OUTPUT VARIABLE(S): ::laieff, ::Pgap_out
1257	!!
1258	!! REFERENCE(S) : Pinty et al; Simplifying the interaction of land surfaces
1259	!! with radiation for relating remote sensing products to climate models
1260	!! JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 111, D02116, doi:10.1029/2005JD005952, 2006
1261	!!
1262	!! Pinty et al; Synergy between 1-D and 3-D radiation transfer models
1263	!! to retrieve vegetation canopy properties from remote sensing data
1264	!! JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 109, D21205, doi:10.1029/2004JD005214, 2004
1265	!!
1266	!! Haverd et al; The Canopy Semi-analytic Pgap and radiative transfer
1267	!! (CanSPART) model: Formulation and application
1268	!! AGRICULTURAL AND FOREST METEOROLOGY, VOL. 160, pp. 14-35, 2012
1269	!!
1270	!! FLOWCHART : None
1271	!! \n
1272	!_ ================================================================================================================================
1273
1274	SUBROUTINE effective_lai(npts, nlevels_loc, z_level_loc, circ_class_biomass, &
1275	circ_class_n, veget_max, sinang, biomass, lai_per_level, laieff, z_array2, &
1276	h_array_out, z_array, Pgap_out, lai_forced)
1277
1278	!! 0 Variable and parameter declaration
1279
1280	!! 0.1 Input variables
1281
1282	INTEGER(i_std),INTENT(IN) :: npts !! Domain size - number of pixels (unitless)
1283	INTEGER(i_std),INTENT(IN) :: nlevels_loc !! The number of physical levels over which we
1284	!! we will calculate the gap probability.
1285	REAL(r_std),DIMENSION(:,:,:),INTENT(IN) :: z_level_loc !! the height of the bottom of each of our levels @tex $(m)$ @endtex
1286	!! note that level 1 is closest to the ground
1287	REAL(r_std), DIMENSION(:,:,:,:,:), &
1288	INTENT(IN) :: circ_class_biomass !! Biomass of the componets of the model
1289	!! tree within a circumference
1290	!! class @tex $(gC ind^{-1})$ @endtex
1291	REAL(r_std), DIMENSION(:,:,:), INTENT(IN) :: circ_class_n !! Number of trees within each circ class
1292	!! class @tex $(m^{-2})$ @endtex
1293	REAL(r_std),DIMENSION(:,:),INTENT(IN) :: veget_max !! Maximum fraction of vegetation type including
1294	!! non-biological fraction (unitless)
1295	REAL(r_std), DIMENSION(:), INTENT(in) :: sinang !! the cosine of the view zenith angle
1296	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! Stand level biomass @tex $(gC m^{-2})$ @endtex
1297	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lai_per_level !! This is the LAI per vertical level
1298	!! @tex $(m^{2} m^{-2})$
1299	REAL(r_std), OPTIONAL :: lai_forced !! ONLY FOR DEBUGGING
1300
1301
1302
1303	!! 0.2 Output variables
1304
1305	REAL(r_std),DIMENSION(:,:,:),INTENT(OUT) :: laieff !! Leaf Area Index Effective converts 3D lai into
1306	!! 1D lai for two stream radiation transfer model
1307	!! @tex $(m^{2} m^{-2})$ @endtex
1308	REAL(r_std),DIMENSION(:,:,:),INTENT(OUT),OPTIONAL &
1309	:: Pgap_out !! The transmission probabilility of light through
1310	!! the canopy. This is per level, not cumulative.
1311	!! (unitless, between 0 and 1)
1312
1313	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc), INTENT(OUT) &
1314	:: z_array !! Height above soil of the Pgap points.
1315	!! @tex $(m)$ @endtex
1316	REAL(r_std),DIMENSION(npts,nvm,nlevels_loc), INTENT(OUT) &
1317	:: z_array2 !! Same as z_array, but one less dimension.
1318	!! @tex $(m)$ @endtex
1319	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc), INTENT(OUT) &
1320	:: h_array_out !! An output of h_array, to use in sechiba
1321
1322	!! 0.3 Modified variables
1323
1324	!! 0.4 Local variables
1325
1326	REAL(r_std), DIMENSION(npts) :: solar_angle !! the view angle of every grid square (radians)
1327	INTEGER(i_std) :: ilevel, ipts, icir !! index (unitless)
1328	INTEGER(i_std) :: ivm, jlevel !! index (unitless)
1329	INTEGER(i_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1330	:: upward_array !! Flag to indicate if the view angle is pointing
1331	!! upwards (1) or downwards (0). (unitless)
1332
1333	! These next values need a bit of caution with the units. Since they are derived
1334	! from circ_class_biomass, it seems they have units per land area, in
1335	! addition to the units given here.
1336	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1337	:: ver_axis_array !! The veritical diameter of the tree canopy.
1338	!! @tex $(m)$ @endtex
1339	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1340	:: hor_axis_array !! The horizontal diameter of the tree canopy.
1341	!! @tex $(m)$ @endtex
1342	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1343	:: h_array !! Height to the top of the tree canopy.
1344	!! @tex $(m)$ @endtex
1345
1346
1347	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1348	:: crown_shadow_h !! The projected shadow of an opaque crown.
1349	!! @tex $(m^2)$ @endtex
1350	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1351	:: trunk_shadow_h !! The projected shadow of the trunk.
1352	!! @tex $(m^2)$ @endtex
1353	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1354	:: crown_area_h !! The area of a horizontal plane cutting
1355	!! through the crown.
1356	!! @tex $(m^2)$ @endtex
1357	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1358	:: sbar_h !! The mean free distance through the crown
1359	!! along the path dictated by the view angle.
1360	!! @tex $(m)$ @endtex
1361	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1362	:: ph_array !! The height of the top of the canopy.
1363	!! @tex $(m)$ @endtex
1364	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1365	:: FAVD_array !! The foliage area volume density of
1366	!! the canopy.
1367	!! @tex $(LAI m^{-3})$ @endtex
1368	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1369	:: DBH_array !! Trunk diameter at breast height.
1370	!! @tex $(m)$ @endtex
1371	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1372	:: G_array !! The mean horizontal projection of
1373	!! unit leaf area (leaf orientation function).
1374	!! @tex $(m^{-2})$ @endtex
1375	! These next two values are not currently used since we set G=0.5
1376	!!$ REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1377	!!$ :: thetaL_mean_array
1378	!!$ REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1379	!!$ :: thetaL_sd_array
1380	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_loc) &
1381	:: Pwc_h !! Crown porosity.
1382	!! (unitless)
1383	REAL(r_std),DIMENSION(npts,nvm,nlevels_loc) :: TrunkShadowL !! Trunk shadow, integrated over the population.
1384	!! @tex $(m^2 m^{-2})$ @endtex
1385	REAL(r_std),DIMENSION(npts,nvm,nlevels_loc) :: CrownAreaL !! The area of a horizontal plane cutting
1386	!! through the crown, averaged over population.
1387	!! @tex $(m^2 m^{-2})$ @endtex
1388	! REAL(r_std),DIMENSION(npts,nvm,nlevels_loc) :: Pbc_trunkL
1389	REAL(r_std),DIMENSION(npts,nvm,nlevels_loc) :: PorousCrownShadowL !! The average shadow cast by a tree taking
1390	!! into account porosity.
1391	!! @tex $(m^2)$ @endtex
1392	REAL(r_std),DIMENSION(npts,nvm,nlevels_loc) :: PgapL !! The probability of finding a gap in the
1393	!! in the canopy.
1394	!! (unitless)
1395	REAL(r_std),DIMENSION(npts,nvm,ncirc) :: CrownVolume_h !! The crown volume of a circ class.
1396	!! @tex $(m^3)$ @endtex
1397	REAL(r_std),DIMENSION(npts,nvm) :: CrownVolumeL !! The average crown volume of the stand.
1398	!! @tex $(m^3)$ @endtex
1399	REAL(r_std),DIMENSION(npts,nvm) :: FAVD !! The foliage area volume density of
1400	!! the canopy.
1401	!! @tex $(LAI m^{-3})$ @endtex
1402	REAL(r_std), DIMENSION(npts,nvm) :: ind_loc !! Density of individuals
1403	!! @tex $(m^{-2})$ @endtex
1404	! The following two variables are ways to speed up the calculation
1405	! by not trying to calculate everything, based on if we have
1406	! biomass or individuals for this PFT on the grid square
1407	LOGICAL, DIMENSION(npts,nvm) :: lcalculate_tree !! For trees
1408	LOGICAL, DIMENSION(npts,nvm) :: lcalculate_nontree !! For non-trees
1409
1410	REAL(r_std),DIMENSION(npts,nvm,ncirc,ndist_types) :: values !! An array that holds canopy properties.
1411	REAL(r_std) :: lai_sum !! The total LAI in the canopy,
1412	!! summed across all levels.
1413	!! @tex $(m^{2} m^{-2})$
1414	REAL(r_std),PARAMETER :: max_laieff=20.0 !! A numerical construct to cap
1415	!! the effective LAI so as to
1416	!! not cause numerical problems.
1417	REAL(r_std) :: min_pgap !! The probability that light makes
1418	!! it through a layer of max_laieff thickness.
1419
1420	!_ ================================================================================================================================
1421
1422	IF (bavard.GE.2) WRITE(numout,*) 'Entering effective_lai'
1423
1424	!! 1. Initialize check for surface area conservation
1425	! Veget_max is a INTENT(in) variable and can therefore
1426	! not be changed during the course of this subroutine
1427	! No need to check whether the subroutine preserves the
1428	! total surface area of the pixel.
1429
1430	!! 2. Initialize
1431	laieff(:,:,:) = zero
1432	PGapL(:,:,:) = un
1433
1434	! we need the total number of individuals
1435	ind_loc(:,:)=SUM(circ_class_n(:,:,:),3)
1436	ind_loc(:,1)=zero
1437
1438
1439	! for each PFT and each grid square, let's determine if we actually need to run
1440	! the calculation
1441	lcalculate_tree(:,:)=.FALSE.
1442	lcalculate_nontree(:,:)=.FALSE.
1443	DO ipts=1,npts
1444	DO ivm=1,nvm
1445	IF(is_tree(ivm))THEN
1446	IF(ind_loc(ipts,ivm) .GT. min_stomate .AND. ind_loc(ipts,ivm) .LT. val_exp*ncirc)THEN
1447	lcalculate_tree(ipts,ivm)=.TRUE.
1448	ENDIF
1449	ELSE
1450	IF(veget_max(ipts,ivm) .GT. min_stomate) lcalculate_nontree(ipts,ivm)=.TRUE.
1451
1452	ENDIF
1453	ENDDO
1454	ENDDO
1455
1456	! for some reason, sinang is really the cosine of the solar zenith
1457	solar_angle(:)=ACOS(sinang(:))
1458
1459	! We need to derive some diagnostic variables from the biomass of each
1460	! tree in the circ distribution, based on allometic relations. We need the
1461	! height of the vegetation, the crown projected surface area, and the trunk
1462	! diameter
1463	CALL derive_biomass_quantities(npts, nvm, ncirc, circ_class_n, &
1464	circ_class_biomass, values)
1465
1466	!************************************
1467	! this is taken directly from Jenny and Vanessa's code
1468	! I've replaced nlayer in their code with nvm, since it seems to me they are
1469	! referring to canopy layers and not strict height levels...I feel these are
1470	! represented by our PFTs, although this is probably better represented by
1471	! species.
1472	! I've also read a little about the WHERE and FORALL constructs, and it seems
1473	! that there is no real advantage to using them over traditional nested loops
1474	! since compilers have gotten really good at loop optimization...I find nested
1475	! loops to be a bit more readable and easier to debug, so I've used those, even
1476	! though the order of the array indices are not necessarily optimized.
1477
1478	DO ipts = 1,npts
1479	DO ivm=1,nvm
1480	IF (.NOT. lcalculate_tree(ipts,ivm)) CYCLE
1481	DO icir=1,ncirc
1482	! we are doing spheres right now, so the horizontal and vertical
1483	! axes have the same distributions
1484	ver_axis_array(ipts,ivm,icir,:) = values(ipts,ivm,icir,icndiaver)
1485	hor_axis_array(ipts,ivm,icir,:) = values(ipts,ivm,icir,icndiahor)
1486	ph_array(ipts,ivm,icir,:) = values(ipts,ivm,icir,iheight)
1487	DBH_array(ipts,ivm,icir,:) = values(ipts,ivm,icir,idiameter)
1488	ENDDO
1489
1490	! these two values are taken from the spherical LAD description in
1491	! section 4.2 of the paper. Looks like they are in degrees.
1492	! We don't need them since we set G=0.5 below
1493	!thetaL_mean_array(ipts,ivm,:,:) = 57.0_r_std
1494	!thetaL_sd_array(ipts,ivm,:,:) = 20.0_r_std
1495
1496	! put some input values into arrays that are used in Vanessa's code
1497
1498	DO icir = 1,ncirc
1499	h_array(ipts,ivm,icir,:) = values(ipts,ivm,icir,iheight)
1500	ENDDO
1501	ENDDO ! loop over PFTS
1502	ENDDO ! end loop over grid
1503
1504	z_array(:,:,:,:) = zero
1505	z_array2(:,:,:) = zero
1506	DO icir=1,ncirc
1507	! these are the z values that we're computing the gap probability for,
1508	! i.e. the bottom of each level. Level 1 is the level next to the soil.
1509	z_array(:,:,icir,:) = z_level_loc(:,:,:)
1510	ENDDO
1511	z_array2(:,:,:) = z_level_loc(:,:,:)
1512
1513	! 0 means we are looking downward, which is what we need
1514	! for the albedo since that is the direction the light is heading
1515	upward_array = 0
1516
1517	DO ipts=1,npts
1518	DO ivm=1,nvm
1519	IF (.NOT. lcalculate_tree(ipts,ivm)) CYCLE
1520	DO icir=1,ncirc
1521	DO ilevel=1,nlevels_loc
1522	! compute the shadow of the crown projected onto each height for
1523	! each view angle using Appendix A in the manuscript...careful,
1524	! there are typos in the manuscript, but these routines came
1525	! directly from Vanessa and Jenny and should be correct if
1526	! either of the elipsoid axes are zero, this routine will crash.
1527	! These limits are arbitrary.
1528	IF(ver_axis_array(ipts,ivm,icir,ilevel) .GT. laieff_zero_cutoff .AND. &
1529	hor_axis_array(ipts,ivm,icir,ilevel) .GT. laieff_zero_cutoff)THEN
1530	crown_shadow_h(ipts,ivm,icir,ilevel) = spheroid_shadow(solar_angle(ipts),&
1531	upward_array(ipts,ivm,icir,ilevel),&
1532	ver_axis_array(ipts,ivm,icir,ilevel),&
1533	hor_axis_array(ipts,ivm,icir,ilevel),&
1534	h_array(ipts,ivm,icir,ilevel),&
1535	z_array(ipts,ivm,icir,ilevel))
1536
1537
1538	! this computes Eq. 10 in the manuscript, the mean path through
1539	! each crown
1540	sbar_h(ipts,ivm,icir,ilevel) = sbar(h_array(ipts,ivm,icir,ilevel),&
1541	ver_axis_array(ipts,ivm,icir,ilevel)/2.0_r_std,&
1542	hor_axis_array(ipts,ivm,icir,ilevel)/2.0_r_std,&
1543	z_array(ipts,ivm,icir,ilevel),&
1544	crown_shadow_h(ipts,ivm,icir,ilevel),&
1545	solar_angle(ipts), upward_array(ipts,ivm,icir,ilevel))
1546	ELSE
1547	crown_shadow_h(ipts,ivm,icir,ilevel)=zero
1548	sbar_h(ipts,ivm,icir,ilevel) = zero
1549	ENDIF ! checking for ver_axis_array(ipts,ivm,icir,ilevel)
1550
1551	! Notice that this crown area is different than the crown area
1552	! we calculate above. This crown area is the area on a plane
1553	! cutting through the crown, whereas that calculated above is
1554	! the crown area projected on the ground from an overhead sun
1555
1556	IF(ver_axis_array(ipts,ivm,icir,ilevel) .GT. laieff_zero_cutoff .AND. &
1557	hor_axis_array(ipts,ivm,icir,ilevel) .GT. laieff_zero_cutoff)THEN
1558	crown_area_h(ipts,ivm,icir,ilevel) = spheroid_area(solar_angle(ipts),&
1559	upward_array(ipts,ivm,icir,ilevel),&
1560	ver_axis_array(ipts,ivm,icir,ilevel),&
1561	hor_axis_array(ipts,ivm,icir,ilevel),&
1562	h_array(ipts,ivm,icir,ilevel),&
1563	z_array(ipts,ivm,icir,ilevel))
1564	ELSE
1565	crown_area_h(ipts,ivm,icir,ilevel)=zero
1566	ENDIF
1567
1568
1569	! This subroutine takes a lot of time, and we're not currently using the trunks
1570	! below. So let's just comment it out for now.
1571	!!$ ! Now compute the shadow of the trunks. We need the same check
1572	!!$ ! as above in case the diameter is zero. No shadow for such a
1573	!!$ ! skinny tree
1574	!!$ IF(DBH_array(ipts,ivm,icir,ilevel) .GT. laieff_zero_cutoff)THEN
1575	!!$ trunk_shadow_h(ipts,ivm,icir,ilevel) = trunk_shadow(solar_angle(ipts),&
1576	!!$ upward_array(ipts,ivm,icir,ilevel),&
1577	!!$ DBH_array(ipts,ivm,icir,ilevel),&
1578	!!$ h_array(ipts,ivm,icir,ilevel),&
1579	!!$ z_array(ipts,ivm,icir,ilevel))
1580	!!$ ELSE
1581	!!$ trunk_shadow_h(ipts,ivm,icir,ilevel)=zero
1582	!!$ ENDIF
1583
1584	ENDDO ! loop over levels
1585	ENDDO ! loop over circ classes
1586
1587	! reassign the crown volume calculated above to an array used in Vanessa's code
1588	DO icir=1,ncirc
1589	CrownVolume_h(ipts,ivm,icir) = values(ipts,ivm,icir,icnvol)
1590	ENDDO
1591
1592	! do some integrations over the height distribution to find average values
1593	! This is equivalent to taking a weighted average by the number of
1594	! individuals in each circ class.
1595
1596	DO ilevel=1,nlevels_loc
1597	!************************
1598	! find the average trunk shadow per square meter
1599	! this is Eq. 5 in the manuscript
1600	! CrownShadowL(ipts,ivm,ilevel)=zero
1601	! in the code, some properties of opaque crowns are calculated, but
1602	! those are not used in the calculation of Pgap so I'm leaving
1603	! them out, like CrownShadowL
1604	!***********************
1605	TrunkShadowL(ipts,ivm,ilevel)=zero
1606	CrownAreaL(ipts,ivm,ilevel)=zero
1607	DO icir=1,ncirc
1608	!++++ CHECK +++++
1609	! Right now we are not using trunks.
1610	! TrunkShadowL(ipts,ivm,ilevel)=TrunkShadowL(ipts,ivm,ilevel)+&
1611	! trunk_shadow_h(ipts,ivm,icir,ilevel)*&
1612	! circ_class_n(ipts,ivm,icir)
1613	!+++++++++++++++++
1614	CrownAreaL(ipts,ivm,ilevel)=CrownAreaL(ipts,ivm,ilevel)+&
1615	crown_area_h(ipts,ivm,icir,ilevel)*&
1616	circ_class_n(ipts,ivm,icir)
1617
1618	ENDDO ! loop over circ classes
1619	CrownAreaL(ipts,ivm,ilevel)=CrownAreaL(ipts,ivm,ilevel)&
1620	/ind_loc(ipts,ivm)
1621
1622	! this is part of Eq. 3 in the manuscript
1623	!++++ CHECK +++++
1624	! Right now we are not using trunks
1625	!TrunkShadowL(ipts,ivm,ilevel)=TrunkShadowL(ipts,ivm,ilevel)*&
1626	! ind_loc(ipts,ivm)
1627	!Pbc_trunkL(ipts,ivm,ilevel)=EXP(-TrunkShadowL(ipts,ivm,ilevel))
1628	!+++++++++++++++
1629
1630	ENDDO! loop over levels
1631
1632	! Compute the average foliage volume density across the height distribution.
1633	! Notice that we do not need the area of the grid square here because we have a
1634	! number density of trees instead of just the number of trees, so grid area
1635	! cancels out in the equation
1636
1637	! We need to integrate over the height distribution. This is equivalent
1638	! to taking a weighted average by the number of individuals in each
1639	! circ class.
1640	CrownVolumeL(ipts,ivm)=zero
1641	DO icir=1,ncirc
1642	CrownVolumeL(ipts,ivm)=CrownVolumeL(ipts,ivm)+&
1643	CrownVolume_h(ipts,ivm,icir)*&
1644	circ_class_n(ipts,ivm,icir)
1645	ENDDO
1646
1647	CrownVolumeL(ipts,ivm)=CrownVolumeL(ipts,ivm)/ind_loc(ipts,ivm)
1648
1649	CALL calculate_favd(ipts,ivm, values(ipts,ivm,:,icnvol), circ_class_biomass, &
1650	!! circ_class_n, biomass, favd, lai=lai_forced)
1651	circ_class_n, biomass, favd)
1652
1653	IF(ld_laieff .AND. ipts == test_grid .AND. ivm == test_pft)THEN
1654	WRITE(numout,*) 'FAVD debugging: ',favd(ipts,ivm)
1655	WRITE(numout,*) ' CrownVolumeL debugging: ',CrownVolumeL(ipts,ivm)
1656	DO icir=1,ncirc
1657	WRITE(numout,'(A,I5,F16.8)') 'Crown diameter: ',icir,values(ipts,ivm,icir,icndiaver)
1658	ENDDO
1659	ENDIF
1660
1661	! generate a larger array of values
1662	FAVD_array(ipts,ivm,:,:) = FAVD(ipts,ivm)
1663
1664	! now we calculate the crown porosity, since the crown is not opaque...Eq. 8
1665	! for this we need the projection of the leaf area in the direction of the beam,
1666	! which is done by the Nilson method (end of Section 4.2 in the manuscript)
1667	! the code also did this with the Ross G function, but we'll take Nilson's.
1668	! the leaf angle distributions (LAD) are currently assumed to be spherical, which
1669	! means the mean leaf inclination is 57Â° and the standard deviation is 20Â°,
1670	! described with a two parameter beta function...this should be explored in more
1671	! depth
1672
1673	DO icir=1,ncirc
1674	DO ilevel=1,nlevels_loc
1675	! G_array(ipts,ivm,icir,ilevel)=&
1676	! NilsonG(thetaL_mean_array(ipts,ivm,icir,ilevel),&
1677	! thetaL_sd_array(ipts,ivm,icir,ilevel),laieff_solar_angle)
1678	! this is done slightly differently from the original code,
1679	! because in the two-stream stream albedo method some coefficients
1680	! are done under the assumption of Ross's G function and a
1681	! spherical leaf distribution, which Pinty et al says makes the
1682	! G function constant and equal to 0.5. Using NilsonG and the
1683	! above parameters gives a value of 0.505, so it's not much
1684	! different but this is good for consistency
1685	G_array(ipts,ivm,icir,ilevel)=0.5_r_std
1686
1687	! I am adding a clumping factor here to account for needles
1688	! clumping to shoots, which increases the porosity of the canopy
1689	! and lets more light through. Notice this is not in the
1690	! original Pgap model.
1691	Pwc_h(ipts,ivm,icir,ilevel)=EXP(-G_array(ipts,ivm,icir,ilevel)*&
1692	FAVD_array(ipts,ivm,icir,ilevel)*sbar_h(ipts,ivm,icir,ilevel)&
1693	*leaf_to_shoot_clumping(ivm))
1694	ENDDO ! loop over levels
1695	ENDDO ! loop over circ classes
1696	ENDDO ! loop over PFTs
1697	ENDDO ! loop over grid points
1698
1699	! Calculate the porous crown shadow integrating over the height distribution.
1700	! This is equivalent
1701	! to taking a weighted average by the number of individuals in each
1702	! circ class.
1703	! Eq. 4, and then calculate Pgap
1704	DO ipts=1,npts
1705	DO ivm=1,nvm
1706
1707	! Here we need to calculate the effective LAI for all land types,
1708	! including grass, crops, and bare soil
1709	IF (is_tree(ivm)) THEN
1710
1711	IF(.NOT. lcalculate_tree(ipts,ivm)) CYCLE
1712
1713	! first for forests
1714	DO ilevel=1,nlevels_loc
1715	PorousCrownShadowL(ipts,ivm,ilevel)=zero
1716	DO icir=1,ncirc
1717	PorousCrownShadowL(ipts,ivm,ilevel)=&
1718	PorousCrownShadowL(ipts,ivm,ilevel)+&
1719	crown_shadow_h(ipts,ivm,icir,ilevel)*&
1720	(1.0_r_std-Pwc_h(ipts,ivm,icir,ilevel))*&
1721	circ_class_n(ipts,ivm,icir)
1722	ENDDO
1723
1724	! These next lines are commented out, but I'm leaving them here
1725	! for clarity. The first line normalises the integration done
1726	! in the previous step so that we have the average porous
1727	! crown shadow over the whole stand. The second line multiplies
1728	! by the density to be used in the exponential of the Poisson
1729	! distribution. It's clear from these two lines that
1730	! they cancel each other out.
1731	!!$ PorousCrownShadowL(ipts,ivm,ilevel)=&
1732	!!$ PorousCrownShadowL(ipts,ivm,ilevel)/ind_loc(ipts,ivm)
1733	!!$ PorousCrownShadowL(ipts,ivm,ilevel)=&
1734	!!$ PorousCrownShadowL(ipts,ivm,ilevel)*ind_loc(ipts,ivm)
1735
1736	! this seems to be Eq. 3 in the manuscript
1737	! PGapL(ipts,ivm,ilevel)=Pbc_trunkL(ipts,ivm,ilevel)*&
1738	! EXP(-PorousCrownShadowL(ipts,ivm,ilevel))
1739	! This line doesn't include the effect of trunks. Since our
1740	! trunks in the albedo calculation aren't distinguishable from
1741	! the leaves, and since we don't want trunks absorbing light, we
1742	! remove the influence of trunks for now.
1743	PGapL(ipts,ivm,ilevel)=&
1744	EXP(-PorousCrownShadowL(ipts,ivm,ilevel))
1745	ENDDO ! loop over levels
1746
1747	ELSEIF(ivm == ibare_sechiba)THEN
1748
1749	! For bare soil the transmission probability is always one. The
1750	! only thing that matters to the albedo is the background reflectance
1751	PGapL(ipts,ivm,:)=un
1752
1753	ELSE
1754
1755	IF(.NOT. lcalculate_nontree(ipts,ivm)) CYCLE
1756
1757	! Now for grass and croplands. We treat these exactly the same
1758	! right now, just a slab of vegetation filling space. For this, we
1759	! need to know the amount of vegetation that is in each level, which
1760	! was calculated in the find_lai_per_level routine. It assumes folliage
1761	! is distributed evenly from the ground to the height of the
1762	! grass/crops. Since the LAI will not be correct for grasslands and
1763	! croplands due to the lack of management in ORCHIDEE, we are
1764	! introducing a tunable parameter to help recover the true albedo values.
1765
1766	DO ilevel=1,nlevels_loc
1767	! This Eq. 1 in the manuscript, assuming G(theta)=0.5 (spherical
1768	! distribution) and using the LAI value with the adjusted LAI
1769	! mentioned above.
1770	! Notice that this lai_temp value has to be cumulative, since that
1771	! of the forests is and we correct Pgap based on that assumption below.
1772	! This means that we need the total LAI through which light has passed
1773	! in order to get here, i.e. the sum of all levels higher than and
1774	! including this level.
1775	!!$ lai_sum = biomass_to_lai(biomass(ipts,ivm,ileaf,icarbon),&
1776	!!$ ivm)
1777	lai_sum=zero
1778	DO jlevel=nlevels_loc,ilevel,-1
1779	lai_sum=lai_sum+lai_per_level(ipts,ivm,jlevel)
1780	ENDDO
1781
1782	PGapL(ipts,ivm,ilevel)=EXP(-0.5_r_stdlai_sum&
1783	lai_correction_factor(ivm)/COS(solar_angle(ipts)))
1784	IF(ipts == test_grid .AND. ivm == test_pft .AND. ld_laieff)THEN
1785	WRITE(numout,*) 'PgapL: ',PGapL(ipts,ivm,ilevel)
1786	WRITE(numout,*) 'lai_sum,lai_correction_factor(ivm): ',lai_sum,lai_correction_factor(ivm)
1787	WRITE(numout,*) 'COS(solar_angle(ipts)),ilevel: ',COS(solar_angle(ipts)),ilevel
1788	ENDIF
1789	ENDDO
1790
1791	ENDIF ! check for tree
1792	ENDDO ! loop over PFTs
1793	ENDDO ! loop over grid points
1794
1795	! Find the pgap for each level seperately...this makes PGap a measure of
1796	! the chance to make it through just this level, as opposed to the
1797	! chance to make it through this level and all the above levels.
1798
1799	! I do not need to CYCLE if lcalculate is not equal to
1800	! true here because I initialize PgapL to be
1801	! equal to 1 on all levels above. The calculation
1802	! of this loop is much less expensive than
1803	! the other loops.
1804	DO ipts=1,npts
1805	DO ivm=1,nvm
1806
1807	IF(ld_laieff)THEN
1808	IF((ipts == test_grid) .AND. (ivm == test_pft))THEN
1809	WRITE(numout,*) 'Cumulative Pgap'
1810	DO ilevel=nlevels_loc,1,-1
1811	WRITE(numout,*) ilevel,PgapL(ipts,ivm,ilevel)
1812	ENDDO
1813	ENDIF
1814	ENDIF
1815
1816	DO ilevel=1,nlevels_loc-1
1817
1818	! include a check for numerical stability.
1819	! The minumum Pgap is a function of the maximum
1820	! effective LAI that we allow. The higher the
1821	! effective LAI, the lower the amount of light
1822	! which will be allowed through the layer. This
1823	! number cannot be so high as to cause numerical
1824	! problems, though. The value of 20 here is
1825	! arbitrary, but doesn't seem to cause any problems.
1826	min_pgap=exp(-max_laieff/2.0/cos(solar_angle(ipts)))
1827
1828	IF(PgapL(ipts,ivm,ilevel+1) .GT. min_pgap )THEN
1829	PgapL(ipts,ivm,ilevel)= PgapL(ipts,ivm,ilevel)/ &
1830	PgapL(ipts,ivm,ilevel+1)
1831	ELSE
1832	! This is a problem. We can't really say anything about
1833	! the gap probability for this level. However, if the
1834	! above level is so dark, no light will make it
1835	! to this level, anyway. So let's just make this
1836	! level opaque. This is equivalent to giving a maximum
1837	! effective LAI value.
1838	PgapL(ipts,ivm,ilevel)=exp(-max_laieff/2.0/cos(solar_angle(ipts)))
1839	ENDIF
1840
1841	! Include a check. PgapL should always be between zero and
1842	! one. Due to numerical issues, it's possible that it won't
1843	! be sometimes. We don't check to see if it's less than
1844	! zero because that will be caught later, and it's a more
1845	! serious issue.
1846	IF(PgapL(ipts,ivm,ilevel) .GT. un) PgapL(ipts,ivm,ilevel)=un
1847
1848	ENDDO ! loop over levels
1849
1850	! compute the effective LAI, based on Eq. 25 in Pintys 2004 paper cited
1851	! in the references above
1852
1853	IF(test_pft == ivm .AND. test_grid == ipts .AND. ld_laieff)THEN
1854	WRITE(numout,*) 'Effective lai'
1855	WRITE(numout,*) 'ipts,ivm,nlevels_loc',ipts,ivm,nlevels_loc
1856	ENDIF
1857
1858	DO ilevel=nlevels_loc,1,-1
1859	! include a check for numerical stability...min_stomate is arbitrary
1860	IF(PgapL(ipts,ivm,ilevel) .GT. min_stomate)THEN
1861	laieff(ilevel,ipts,ivm)=-2.0_r_stdCOS(solar_angle(ipts))&
1862	LOG(PgapL(ipts,ivm,ilevel))
1863	ELSE
1864	laieff(ilevel,ipts,ivm)=-2.0_r_stdCOS(solar_angle(ipts))&
1865	LOG(min_stomate)
1866	ENDIF
1867
1868	! If the value is small enough or negative, we set it equal to zero.
1869	! A negative value for the effective LAI does not make sense in the
1870	! context of Eq. 12 in Pinty 2006. If the effective LAI is zero, the
1871	! uncollided transmission will be one. If LAI is positive, the
1872	! uncollided transmission will be between zero and one, since light
1873	! is intercepted by the canopy elements. A negative effective
1874	! LAI results in uncollided transmitted light being greater than one,
1875	! which means there is a light source in the canopy.
1876	IF(laieff(ilevel,ipts,ivm) .LE. min_stomate)THEN
1877	IF( laieff(ilevel,ipts,ivm) .LE. -min_stomate) THEN
1878	WRITE(numout,*) 'ilevel,ipts,ivm', ilevel,ipts,ivm
1879	WRITE(numout,*) 'laieff(ilevel,ipts,ivm)', laieff(ilevel,ipts,ivm)
1880	CALL ipslerr_p (3,'effective_lai', 'We should not have a negative effective LAI.','','')
1881	ENDIF
1882	laieff(ilevel,ipts,ivm) = zero
1883	ENDIF
1884
1885
1886	ENDDO ! loop over levels
1887
1888	IF(test_pft == ivm .AND. test_grid == ipts .AND. ld_laieff)THEN
1889	WRITE(numout,*) 'ilevel, laieff, PgapL: '
1890	DO ilevel=nlevels_loc,1,-1
1891	WRITE(numout,'(I5,2E20.10)') ilevel,laieff(ilevel,ipts,ivm),PgapL(ipts,ivm,ilevel)
1892	ENDDO
1893	ENDIF
1894
1895	ENDDO ! loop over PFTs
1896	ENDDO ! loop over grid points
1897
1898
1899	h_array_out = h_array
1900
1901	IF(PRESENT(Pgap_out))THEN
1902	Pgap_out(:,:,:)=PgapL(:,:,:)
1903	ENDIF
1904
1905	IF (bavard.GE.2) WRITE(numout,*) 'Leaving effective_lai'
1906
1907	END SUBROUTINE effective_lai
1908
1909
1910
1911	!! ================================================================================================================================
1912	!! FUNCTION : spheroid_shadow
1913	!!
1914	!>\BRIEF Calculates the projected shadow of a tree canopy at different heights and
1915	!! sun angles, assuming a spheriod canopy
1916	!!
1917	!! DESCRIPTION : Taken with permission from V. HAVERD's code, based on the equations in Appendix A
1918	!! of her paper
1919	!!
1920	!! RECENT CHANGE(S) : None
1921	!!
1922	!! MAIN OUTPUT VARIABLE(S): ::spheroid_shadow
1923	!!
1924	!! REFERENCE(S) :
1925	!! Haverd et al; The Canopy Semi-analytic Pgap and radiative transfer
1926	!! (CanSPART) model: Formulation and application, AGRICULTURAL AND FOREST METEOROLOGY
1927	!! Vol. 160, pp. 14-35, 2012
1928	!!
1929	!! FLOWCHART : None
1930	!! \n
1931	!_ ================================================================================================================================
1932
1933	REAL(r_std) FUNCTION spheroid_shadow(theta,upward,T,D,&
1934	H,z)
1935
1936	!! 0 Variable and parameter declaration
1937
1938	!! 0.1 Input variables
1939	REAL(r_std), INTENT(IN) :: theta ! view angle (radians)
1940	REAL(r_std), INTENT(IN) :: T ! vertical diameter of the spheroid (m)
1941	REAL(r_std), INTENT(IN) :: D ! horizontal diameter (m)
1942	REAL(r_std), INTENT(IN) :: H ! tree height (top of crown)
1943	REAL(r_std), INTENT(IN) :: z ! height above ground
1944	INTEGER(i_std), INTENT(IN) :: upward ! ==1 for upward view, ==0 for downward view
1945
1946
1947	!! 0.2 Output variables
1948	! REAL(r_std) :: spheroid_shadow
1949
1950
1951	!! 0.3 Modified variables
1952
1953	!! 0.4 Local variables
1954	REAL(r_std) :: thetapr, Ypr, Y1pr, Y1, X1, D1pr, D2pr, Y
1955
1956	!_ ================================================================================================================================
1957
1958	thetapr = ATAN(T/D*TAN(theta)) ! zenith angle in spherical space
1959
1960
1961	! the conventions are described in Appendix A of this paper
1962	IF (upward==1) THEN
1963	Y = z-H+T/2
1964	ELSE
1965	Y = H-T/2-z
1966	ENDIF
1967	Ypr = Y*D/T
1968	Y1pr = D/2.0TAN(thetapr)/SQRT(1+TAN(thetapr)*2)
1969	Y1 = T/D*Y1pr
1970	IF (TAN(thetapr).NE.0.0) THEN
1971	X1 = Y1pr/TAN(thetapr)
1972	ELSE
1973	X1 = D/2
1974	END IF
1975
1976
1977
1978	! first calculate major axis (D1pr) of elliptic shadow
1979
1980	IF (Y<=-T/2.0) THEN ! below crown bottom
1981	D1pr = 0.0
1982	END IF
1983
1984	IF (Y>-T/2 .AND. Y<=-Y1) THEN !between crown bottom and lower tangent point
1985	D1pr = DSQRT(1.0-Ypr24.0/D**2)
1986	END IF
1987
1988	IF (Y>-Y1 .AND. Y <Y1) THEN !between lower and upper tangent points
1989	D1pr = D/2.0SQRT(1.0-MIN(Ypr24.0/D**2,1.0)) +X1 - &
1990	TAN(thetapr)*(-Y1pr-Ypr)
1991	END IF
1992
1993	IF (Y > Y1) THEN !above upper tangent point
1994	D1pr = D/COS(thetapr)
1995	END IF
1996
1997	! second calculate minor axis (D2pr) of elliptic shadow
1998	! (this is just the diameter of the circular cross-section at height z)
1999	IF (Y <= -T/2) THEN
2000	D2pr = 0.0 !shadow breadth
2001	END IF
2002
2003	IF (Y<=0 .AND. Y >-T/2) THEN
2004	D2pr = 2.0SQRT((D/2.0)2-Ypr*2) !shadow breadth
2005	END IF
2006
2007	IF (Y>=0) THEN
2008	D2pr = D !shadow breadth
2009	END IF
2010
2011
2012	! shadow area
2013	spheroid_shadow = pi/4.0D2prD1pr
2014
2015	END FUNCTION spheroid_shadow
2016
2017
2018	!! ================================================================================================================================
2019	!! FUNCTION : trunk_shadow
2020	!!
2021	!>\BRIEF ! Calculates the projected shadow of a tree trunk at different heights and
2022	!! sun angles, assuming a cone shape
2023	!!
2024	!! DESCRIPTION : Taken with permission from V. HAVERD's code, based on the equations in Appendix E
2025	!! of her paper
2026	!!
2027	!! RECENT CHANGE(S): None
2028	!!
2029	!! RETURN VALUE : trunk_shadow
2030	!!
2031	!! REFERENCE(S) : Haverd et al; The Canopy Semi-analytic Pgap and radiative transfer
2032	!! (CanSPART) model: Formulation and application, AGRICULTURAL AND FOREST METEOROLOGY
2033	!! Vol. 160, pp. 14-35, 2012
2034	!!
2035	!! FLOWCHART : None
2036	!! \n
2037	!_ ================================================================================================================================
2038
2039	REAL(r_std) FUNCTION trunk_shadow(theta,upward,D,H,z)
2040	!! 0 Variable and parameter declaration
2041
2042	!! 0.1 Input variables
2043	REAL(r_std), INTENT(IN) :: theta ! view angle @tex $(radians)$ @endtex
2044	REAL(r_std), INTENT(IN) :: D ! trunk diameter @tex $(m)$ @endtex
2045	REAL(r_std), INTENT(IN) :: H ! tree height (top of crown) @tex $(m)$ @endtex
2046	REAL(r_std), INTENT(IN) :: z ! height above ground @tex $(m)$ @endtex
2047	INTEGER(i_std), INTENT(IN) :: upward ! ==1 for upward view, ==0 for downward view
2048
2049	!! 0.2 Output variables
2050
2051	!! 0.3 Modified variables
2052
2053	!! 0.4 Local variables
2054	REAL(r_std) :: sinalpha, alpha
2055
2056	!_ ================================================================================================================================
2057
2058	sinalpha = D/(2.HTAN(theta))
2059	IF (ABS(sinalpha).LT.1.0.AND.D.GT.0.0) THEN
2060	alpha = 2.*ASIN(sinalpha)
2061
2062	IF (upward==1) THEN
2063	! shadow is bounded by base-of-trunk circle; trunk cross-section at z
2064	! and tangents joining the two circles
2065	IF (z<H) THEN
2066	trunk_shadow = 1./4. * D*2(1.-(1.-z/H)**2)/(TAN(alpha/2.)) + &
2067	1./8.D2(pi+alpha) + 1./8.D2(1-z/H)*2(pi-alpha)
2068	ELSE
2069	trunk_shadow = 1./4. * D*21./(TAN(alpha/2.)) + 1./8.D2(pi+alpha)
2070	ENDIF
2071
2072	trunk_shadow = trunk_shadow - pi * (D**2)/4.
2073	ELSEIF (upward==0) THEN
2074
2075	! shadow is bounded by trunk cross-section at z and tangents joining two circles
2076	IF (z<H) THEN
2077	trunk_shadow = 1./4. * D*2(1.-z/H)*21./(TAN(alpha/2.)) + &
2078	1./8.D2(1.-z/H)*2(pi+alpha)
2079	ELSE
2080	trunk_shadow = 0.
2081	ENDIF
2082
2083	ENDIF
2084
2085	ELSE
2086	trunk_shadow = 0.0
2087	ENDIF
2088
2089	END FUNCTION trunk_shadow
2090
2091
2092	!! ================================================================================================================================
2093	!! FUNCTION : sbar
2094	!!
2095	!>\BRIEF ! Calculates the mean free distance through the crown, Eq. 10 in Haverd et al.
2096	!!
2097	!! DESCRIPTION : Taken with permission from V. HAVERD's code
2098	!!
2099	!! RECENT CHANGE(S): None
2100	!!
2101	!! RETURN VALUE : sbar
2102	!!
2103	!! REFERENCE(S) : Haverd et al; The Canopy Semi-analytic Pgap and radiative transfer
2104	!! (CanSPART) model: Formulation and application, AGRICULTURAL AND FOREST METEOROLOGY
2105	!! Vol. 160, pp. 14-35, 2012
2106	!!
2107	!! FLOWCHART : None
2108	!! \n
2109	!_ ================================================================================================================================
2110
2111	REAL(r_std) FUNCTION sbar(x, b, r, z, area,&
2112	theta, upward )
2113	!! 0 Variable and parameter declaration
2114
2115	!! 0.1 Input variables
2116	REAL(r_std), INTENT(IN) :: x ! tree height
2117	! @tex $(m)$ @endtex
2118	REAL(r_std), INTENT(IN) :: b ! vertical crown radius
2119	! @tex $(m)$ @endtex
2120	REAL(r_std), INTENT(IN) :: r ! horizontal crown radius
2121	! @tex $(m)$ @endtex
2122	REAL(r_std), INTENT(IN) :: z ! height above ground
2123	! @tex $(m)$ @endtex
2124	REAL(r_std), INTENT(IN) :: theta ! view angle
2125	! @tex $(radians)$ @endtex
2126	REAL(r_std), INTENT(IN) :: area ! shadow area of partial spheroid
2127	! @tex $(m^2)$ @endtex
2128	INTEGER(i_std), INTENT(IN) :: upward ! ==1 for upward view, ==0
2129	! for downward view
2130
2131	!! 0.2 Output variables
2132
2133	!! 0.3 Modified variables
2134
2135	!! 0.4 Local variables
2136	REAL(r_std) :: Hpr, Y, V_ps
2137
2138	!_ ================================================================================================================================
2139
2140	Hpr = x-b ! Hpr is height at crown centre
2141	IF (upward==1) THEN
2142	Y = z-Hpr ! Y=0 corresponds to crown centre
2143	ELSE
2144	Y = Hpr - z
2145	ENDIF
2146
2147	! initialize the output variables, just in case there is a numerical reason
2148	! none of the loops is hit
2149	V_ps =zero
2150	sbar =zero
2151
2152	IF (Y<b.AND.Y>-b) THEN
2153
2154	! volume of partial spheroid below height z
2155	V_ps = 4.0_r_std/3.0_r_stdpi(r)*2(b) - ((2.0_r_std* b)/3.0_r_std - Y &
2156	+ Y*3/(3.0_r_std b*2)) pi* r**2
2157
2158	! there are some numerical issues here occasionally that give a negative
2159	! extremely small value, which in turn gives a negative sbar which causes
2160	! problems later
2161	IF(ABS(V_ps) .LT. 0.0000001_r_std) V_ps= zero
2162	sbar = V_ps/COS(theta)/area
2163
2164	ENDIF
2165
2166	!!!!!! had to add the check for area
2167	IF (Y>=b .AND. area .NE. zero) THEN
2168	V_ps = 4.0_r_std/3.0_r_stdpi(r)*2(b) ; ! volume of whole spheroid
2169	sbar= V_ps/COS(theta)/area
2170	ENDIF
2171
2172
2173	END FUNCTION sbar
2174
2175	!! ================================================================================================================================
2176	!! FUNCTION : spheroid_area
2177	!!
2178	!>\BRIEF ! Calculates the area of the cross section of a spheroid at a height of z
2179	!!
2180	!! DESCRIPTION : Taken with permission from V. HAVERD's code
2181	!!
2182	!! RECENT CHANGE(S): None
2183	!!
2184	!! RETURN VALUE : spheroid_area
2185	!!
2186	!! REFERENCE(S) : Haverd et al; The Canopy Semi-analytic Pgap and radiative transfer
2187	!! (CanSPART) model: Formulation and application, AGRICULTURAL AND FOREST METEOROLOGY
2188	!! Vol. 160, pp. 14-35, 2012
2189	!!
2190	!! FLOWCHART : None
2191	!! \n
2192	!_ ================================================================================================================================
2193
2194	REAL(r_std) FUNCTION spheroid_area(theta,upward,T,D,H,z)
2195
2196	!! 0 Variable and parameter declaration
2197
2198	!! 0.1 Input variables
2199	REAL(r_std), INTENT(IN) :: theta ! view angle
2200	! @tex $(radians)$ @endtex
2201	REAL(r_std), INTENT(IN) :: T ! vertical diameter
2202	! @tex $(m)$ @endtex
2203	REAL(r_std), INTENT(IN) :: D ! horizontal diameter
2204	! @tex $(m)$ @endtex
2205	REAL(r_std), INTENT(IN) :: H ! tree height (top of crown)
2206	! @tex $(m)$ @endtex
2207	REAL(r_std), INTENT(IN) :: z ! height above ground
2208	! @tex $(m)$ @endtex
2209	INTEGER(i_std), INTENT(IN) :: upward ! ==1 for upward view,
2210	! ==0 for downward view
2211
2212
2213	!! 0.2 Output variables
2214
2215	!! 0.3 Modified variables
2216
2217	!! 0.4 Local variables
2218	REAL(r_std) :: thetapr, Ypr, Y1pr, Y1, X1, D2pr, Y
2219
2220	!_ ================================================================================================================================
2221
2222	thetapr = ATAN(T/D*TAN(theta)) ! zenith angle in spherical space
2223
2224
2225	IF (upward==1) THEN
2226	Y = z-H+T/2
2227	ELSE
2228	Y = H-T/2-z
2229	ENDIF
2230	Ypr = Y*D/T
2231	Y1pr = D/2.0TAN(thetapr)/SQRT(1+TAN(thetapr)*2)
2232	Y1 = T/D*Y1pr
2233	IF (TAN(thetapr).NE.0.0) THEN
2234	X1 = Y1pr/TAN(thetapr)
2235	ELSE
2236	X1 = D/2
2237	END IF
2238
2239
2240	! calculate minor axis (D2pr) of elliptic shadow
2241	! (this is just the diameter of the circular cross-section at height z)
2242
2243
2244	IF (Y<=T/2 .AND. Y >-T/2) THEN
2245	D2pr = 2.0SQRT((D/2.0)2-Ypr*2) !shadow breadth
2246	ELSE
2247	D2pr = 0 !shadow breadth
2248	END IF
2249
2250
2251	! shadow area
2252	spheroid_area = pi/4.0D2pr*2
2253
2254	END FUNCTION spheroid_area
2255
2256	!! ================================================================================================================================
2257	!! FUNCTION : NilsonG
2258	!!
2259	!>\BRIEF ! Calculates the projection of unit foliage area G(theta) as described by
2260	!! Wang et al, 2007, referenced in section 4.2 from V. Haverd's paper
2261	!!
2262	!! DESCRIPTION : Taken with permission from V. HAVERD's code
2263	!!
2264	!! RECENT CHANGE(S): None
2265	!!
2266	!! RETURN VALUE : NilsonG
2267	!!
2268	!! REFERENCE(S) : Haverd et al; The Canopy Semi-analytic Pgap and radiative transfer
2269	!! (CanSPART) model: Formulation and application, AGRICULTURAL AND FOREST METEOROLOGY
2270	!! Vol. 160, pp. 14-35, 2012
2271	!!
2272	!! FLOWCHART : None
2273	!! \n
2274	!_ ================================================================================================================================
2275
2276	REAL(r_std) FUNCTION NilsonG(mean_thetaL,sd_thetaL,theta_deg)
2277	!! 0 Variable and parameter declaration
2278
2279	!! 0.1 Input variables
2280	REAL(r_std), INTENT(IN) :: mean_thetaL ! the mean leaf angle inclination
2281	! @tex $(degrees)$ @endtex
2282	REAL(r_std), INTENT(IN) :: sd_thetaL ! the standard deviation of the leaf angle inclination
2283	! @tex $(degrees)$ @endtex
2284	REAL(r_std), INTENT(IN) :: theta_deg ! the solar zenith angle
2285	! @tex $(degrees)$ @endtex
2286
2287	!! 0.2 Output variables
2288
2289	!! 0.3 Modified variables
2290
2291	!! 0.4 Local variables
2292	REAL(r_std) :: var_t, thl_bar, theta, thetaL, psi
2293	REAL(r_std), ALLOCATABLE :: A(:), f(:)
2294	INTEGER(i_std) :: k, nthetaL
2295	!_ ================================================================================================================================
2296
2297	var_t = (sd_thetaL/90.0_r_std)**2
2298	thl_bar = mean_thetaL*pi/180.0_r_std
2299	theta = theta_deg*pi/180.0_r_std
2300	nthetaL = 20
2301	ALLOCATE (A(nthetaL))
2302	ALLOCATE (f(nthetaL))
2303	DO k=1,nthetaL
2304	thetaL = 0.0_r_std + 89.0_r_std/nthetaL(k)pi/180.0_r_std
2305	IF (ABS((1/TAN(theta))*(1/TAN(thetaL)))>1) THEN
2306	A(k) = COS(theta)*COS(thetaL);
2307	ELSE
2308	psi = ACOS((1.0_r_std/TAN(theta))*(1.0_r_std/TAN(thetaL)));
2309	A(k) = COS(theta)COS(thetaL)(1.0_r_std+2.0_r_std/pi*(TAN(psi)-psi));
2310	ENDIF
2311	f(k) = BETA_LAD(thl_bar,var_t, thetaL)
2312	ENDDO
2313	NilsonG = SUM(A(1:nthetaL)*f(1:nthetaL))/SUM(f(1:nthetaL))
2314
2315	DEALLOCATE(A)
2316	DEALLOCATE(F)
2317
2318	END FUNCTION NilsonG
2319
2320	!! ================================================================================================================================
2321	!! FUNCTION : BETA_LAD
2322	!!
2323	!>\BRIEF ! Calculates the beta leaf angle distribution as described by
2324	!! Wang et al, 2007, referenced in section 4.2 from V. Haverd's paper
2325	!!
2326	!! DESCRIPTION : Taken with permission from V. HAVERD's code
2327	!!
2328	!! RECENT CHANGE(S): None
2329	!!
2330	!! RETURN VALUE : BETA_LAD
2331	!!
2332	!! REFERENCE(S) : Haverd et al; The Canopy Semi-analytic Pgap and radiative transfer
2333	!! (CanSPART) model: Formulation and application, AGRICULTURAL AND FOREST METEOROLOGY
2334	!! Vol. 160, pp. 14-35, 2012
2335	!!
2336	!! FLOWCHART : None
2337	!! \n
2338	!_ ================================================================================================================================
2339
2340	REAL(r_std) FUNCTION BETA_LAD(thl_bar,var_t, theta)
2341
2342	!! 0 Variable and parameter declaration
2343
2344	!! 0.1 Input variables
2345	REAL(r_std), INTENT(IN):: thl_bar,var_t, theta
2346
2347	!! 0.2 Output variables
2348
2349	!! 0.3 Modified variables
2350
2351	!! 0.4 Local variables
2352	REAL(r_std) :: tbar, sig0_2, sigt_2, mu, nu, t, num
2353	!_ ================================================================================================================================
2354
2355	tbar = 2.0_r_std*thl_bar/(pi)
2356	sig0_2 = tbar*(1.0_r_std-tbar)
2357	sigt_2 = var_t
2358	num = sig0_2/sigt_2 - 1.0_r_std
2359	! Test for negative num as this causes distribution to fail.
2360	IF (num <= zero) THEN
2361	num = 0.01_r_std
2362	!write(,)'WARNING: leaf angle SD too large'
2363	ENDIF
2364
2365	nu = tbar*num
2366	mu = (1.0_r_std-tbar)*num
2367
2368	t = 2.0_r_std*theta/pi
2369	beta_lad = ((1.0_r_std-t)*(mu-1.0_r_std))(t**(nu-1.0_r_std))/beta(mu,nu)
2370
2371	END FUNCTION BETA_LAD
2372	!********************************************************************************
2373
2374	!! ================================================================================================================================
2375	!! FUNCTION : beta
2376	!!
2377	!>\BRIEF ! Calculates the beta leaf angle distribution as described by
2378	!! Wang et al, 2007, referenced in section 4.2 from V. Haverd's paper
2379	!!
2380	!! DESCRIPTION : Taken with permission from V. HAVERD's code
2381	!!
2382	!! RECENT CHANGE(S): None
2383	!!
2384	!! RETURN VALUE : BETA
2385	!!
2386	!! REFERENCE(S) : Haverd et al; The Canopy Semi-analytic Pgap and radiative transfer
2387	!! (CanSPART) model: Formulation and application, AGRICULTURAL AND FOREST METEOROLOGY
2388	!! Vol. 160, pp. 14-35, 2012
2389	!!
2390	!! FLOWCHART : None
2391	!! \n
2392	!_ ================================================================================================================================
2393
2394	REAL(r_std) FUNCTION beta(z,w)
2395	!! 0 Variable and parameter declaration
2396
2397	!! 0.1 Input variables
2398	REAL(r_std), INTENT(IN):: z,w
2399
2400	!! 0.2 Output variables
2401
2402	!! 0.3 Modified variables
2403
2404	!! 0.4 Local variables
2405	!_ ================================================================================================================================
2406
2407	beta=EXP(gammln(z)+gammln(w)-gammln(z+w))
2408
2409	END FUNCTION beta
2410	!! ================================================================================================================================
2411	!! FUNCTION : gammaln
2412	!!
2413	!>\BRIEF ! Calculates the beta leaf angle distribution as described by
2414	!! Wang et al, 2007, referenced in section 4.2 from V. Haverd's paper
2415	!!
2416	!! DESCRIPTION : Taken with permission from V. HAVERD's code
2417	!!
2418	!! RECENT CHANGE(S): None
2419	!!
2420	!! RETURN VALUE : gammaln
2421	!!
2422	!! REFERENCE(S) : Haverd et al; The Canopy Semi-analytic Pgap and radiative transfer
2423	!! (CanSPART) model: Formulation and application, AGRICULTURAL AND FOREST METEOROLOGY
2424	!! Vol. 160, pp. 14-35, 2012
2425	!!
2426	!! FLOWCHART : None
2427	!! \n
2428	!_ ================================================================================================================================
2429
2430	REAL(r_std) FUNCTION gammln(xx)
2431	!! 0 Variable and parameter declaration
2432
2433	!! 0.1 Input variables
2434	REAL(r_std), INTENT(IN):: xx
2435
2436	!! 0.2 Output variables
2437
2438	!! 0.3 Modified variables
2439
2440	!! 0.4 Local variables
2441	REAL(r_std) :: tmp,x,stp,temp,first,increment
2442	INTEGER(i_std) :: k,k2,n
2443	INTEGER(i_std), PARAMETER :: NPAR_ARTH=16,NPAR2_ARTH=8
2444	REAL(r_std),DIMENSION(6) :: coef,ARTH
2445	!_ ================================================================================================================================
2446
2447	stp = 2.5066282746310005_r_std
2448	coef = (/76.18009172947146_r_std,&
2449	-86.50532032941677_r_std,24.01409824083091_r_std,&
2450	-1.231739572450155_r_std,0.1208650973866179e-2_r_std,&
2451	-0.5395239384953e-5_r_std/)
2452	x=xx
2453	tmp=x+5.5_r_std
2454	tmp=(x+0.5_r_std)*LOG(tmp)-tmp
2455
2456	n = SIZE(coef)
2457	increment = 1.0_r_std
2458	first = x+1.0_r_std
2459	IF (n > 0) arth(1)=first
2460	IF (n <= NPAR_ARTH) THEN
2461	DO k=2,n
2462	ARTH(k)=ARTH(k-1)+increment
2463	END DO
2464	ELSE
2465	DO k=2,NPAR2_ARTH
2466	ARTH(k)=ARTH(k-1)+increment
2467	END DO
2468	temp=increment*NPAR2_ARTH
2469	k=NPAR2_ARTH
2470	DO
2471	IF (k >= n) EXIT
2472	k2=k+k
2473	arth(k+1:MIN(k2,n))=temp+arth(1:MIN(k,n-k))
2474	temp=temp+temp
2475	k=k2
2476	END DO
2477	END IF
2478
2479	gammln=tmp+LOG(stp*(1.000000000190015_r_std+&
2480	SUM(coef(:)/arth(:)))/x)
2481	END FUNCTION gammln
2482
2483	!! ================================================================================================================================
2484	!! FUNCTION : partial_spheroid_vol
2485	!!
2486	!>\BRIEF ! Calculates the beta leaf angle distribution as described by
2487	!! Wang et al, 2007, referenced in section 4.2 from V. Haverd's paper
2488	!!
2489	!! DESCRIPTION : Taken with permission from V. HAVERD's code
2490	!!
2491	!! RECENT CHANGE(S): None
2492	!!
2493	!! RETURN VALUE : partial_spheroid_vol
2494	!!
2495	!! REFERENCE(S) : Haverd et al; The Canopy Semi-analytic Pgap and radiative transfer
2496	!! (CanSPART) model: Formulation and application, AGRICULTURAL AND FOREST METEOROLOGY
2497	!! Vol. 160, pp. 14-35, 2012
2498	!!
2499	!! FLOWCHART : None
2500	!! \n
2501	!_ ================================================================================================================================
2502
2503	REAL(r_std) FUNCTION partial_spheroid_vol(Z_up, Z_down, b, r, h)
2504	!! 0 Variable and parameter declaration
2505
2506	!! 0.1 Input variables
2507	! The units of all these variables can be anything, so long as they
2508	! are all the same.
2509	! I assume
2510	REAL(r_std), INTENT(IN) :: Z_up ! The height of the top cutting plane
2511	REAL(r_std), INTENT(IN) :: Z_down ! The height of the lower cutting plane
2512	REAL(r_std), INTENT(IN) :: b ! Distance from center to top of the spheroid
2513	REAL(r_std), INTENT(IN) :: r ! Distance from center to outside on the XY plane
2514	REAL(r_std), INTENT(IN) :: h ! Height of the top of the spheroid above the ground
2515	!! 0.2 Output variables
2516
2517	!! 0.3 Modified variables
2518
2519	!! 0.4 Local variables
2520	REAL(r_std) :: Z1, Z2, center
2521
2522	!_ ================================================================================================================================
2523
2524	! initialize the output variables
2525	partial_spheroid_vol =zero
2526
2527	! It's possible that none of this spheroid is in this level!
2528	center = h-b
2529
2530	IF(Z_down - center .GE. b)RETURN
2531	IF(Z_up - center .LE. -b)RETURN
2532
2533	! Z1 and Z2 have to be within the limits of the spheroid
2534	! I also want the center at the origin.
2535	Z2=MIN(Z_up-center,b)
2536	Z1=MAX(Z_down-center,-b)
2537
2538	! With a little calculus, I computed the following formula, which
2539	! gives the correct result in the case of the full spheroid,
2540	! i.e. Z2=b and Z1=-b, as well as for a full sphere (r=b)
2541	partial_spheroid_vol = pir2(Z2-Z1-(Z23-Z13)/(3.0_r_stdb*2))
2542
2543	! there are some numerical issues here occasionally that give a negative
2544	! extremely small value, which in turn gives a negative partial_spheroid_vol which causes
2545	! problems later
2546	IF(ABS(partial_spheroid_vol) .LT. min_stomate) partial_spheroid_vol= zero
2547
2548	END FUNCTION partial_spheroid_vol
2549
2550	!! ================================================================================================================================
2551	!! SUBROUTINE : calculate_favd
2552	!!
2553	!>\BRIEF Calculates the foliage area volume density of a canopy.
2554	!!
2555	!! DESCRIPTION : WARNING: This should only be used for trees. Grasses
2556	!! and crops have no canopy volume currently in the model.
2557	!!
2558	!! NOTE:
2559	!!
2560	!! RECENT CHANGE(S) : None
2561	!!
2562	!! MAIN OUTPUT VARIABLE(S): ::favd
2563	!!
2564	!! REFERENCE(S) :
2565	!!
2566	!! FLOWCHART : None
2567	!! \n
2568	!_ ================================================================================================================================
2569
2570	SUBROUTINE calculate_favd(ipts,ivm, cn_vol, circ_class_biomass, &
2571	circ_class_n, biomass, favd)
2572
2573	!! 0 Variable and parameter declaration
2574
2575	!! 0.1 Input variables
2576
2577	INTEGER(i_std),INTENT(IN) :: ipts
2578	INTEGER(i_std),INTENT(IN) :: ivm
2579	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(IN) :: circ_class_biomass !! Biomass of the components of the model
2580	!! tree within a circumference class
2581	!! @tex $(gC ind^{-1})$ @endtex
2582	REAL(r_std), DIMENSION(:,:,:), INTENT(IN) :: circ_class_n !! Number of trees within each circ class
2583	!! @tex $(ind m^{-2})$ @endtex
2584	REAL(r_std), DIMENSION(:), INTENT(IN) :: cn_vol !! Canopy volume for each circ class
2585	!! @tex $(m^{-3} ind^{-1})$ @endtex
2586	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! Stand level biomass
2587	!! @tex $(gC m^{-2})$ @endtex
2588
2589	!! 0.2 Output variables
2590	REAL(r_std), DIMENSION(:,:), INTENT(OUT) :: favd !! The LAI per canopy volume density
2591	!! @tex $(m^{-3} m^{-2})$ @endtex
2592
2593	!! 0.3 Modified variables
2594
2595	!! 0.4 Local variables
2596	INTEGER :: icir
2597	REAL(r_std) :: total_bm
2598	REAL(r_std) :: total_vol
2599	REAL(r_std) :: total_la
2600	REAL(r_std) :: max_height
2601	REAL(r_std) :: lai_temp
2602	REAL(r_std) :: ind_loc
2603	!_ ================================================================================================================================
2604
2605	IF (bavard.GE.2) WRITE(numout,*) 'Entering calculate_favd'
2606
2607	ind_loc = SUM(circ_class_n(ipts,ivm,:))
2608
2609	IF(is_tree(ivm))THEN
2610	! How much total leaf biomass do we have?
2611	! How much total canopy volume do we have?
2612	total_bm=zero
2613	total_vol=zero
2614	! Notice that we don't have to take the surface area of the stand
2615	! into account. total_bm will be in units of gC / m**2, while
2616	! the canopy volume should be in units of m3 / m2. The surface
2617	! area would cancel out when we take the ratio.
2618	DO icir=1,ncirc
2619	total_bm=total_bm+circ_class_biomass(ipts,ivm,icir,ileaf,icarbon)*&
2620	! circ_class_n(ipts,ivm,icir)*surface_area
2621	circ_class_n(ipts,ivm,icir)
2622	total_vol=total_vol+cn_vol(icir)*&
2623	circ_class_n(ipts,ivm,icir)
2624	! circ_class_n(ipts,ivm,icir)*surface_area
2625	IF(ipts == test_grid .AND. ivm == test_pft .AND. ld_laieff)THEN
2626	WRITE(numout,'(A,I5,3F16.8)') 'icir,cn_vol,circ_class_biomass,circ_class_n: ',&
2627	icir,cn_vol(icir),circ_class_biomass(ipts,ivm,icir,ileaf,icarbon),circ_class_n(ipts,ivm,icir)
2628	ENDIF
2629	ENDDO
2630
2631	! Right now, total_vol is in m*3/tree tree/m**2
2632
2633	! How much leaf area do we have on our grid square?
2634	total_la=total_bm*sla(ivm)
2635
2636	! Total_la is unitless. m2 / m2.
2637
2638	IF(ipts == test_grid .AND. ivm == test_pft .AND. ld_laieff)THEN
2639	WRITE(numout,'(A,3F16.8)') &
2640	'total_la,total_vol: ',total_la,total_vol,total_la/total_vol
2641	ENDIF
2642
2643	! What is the foliage area volume density?
2644	IF (total_vol .GT. min_stomate) THEN
2645	! 1 / m3 (canopy volume) / m2 (land area)
2646	FAVD(ipts,ivm)=total_la/total_vol
2647	ELSE
2648	FAVD(ipts,ivm)=zero
2649	ENDIF
2650
2651
2652	ELSE
2653
2654	! What is the height of our grass/crop?
2655
2656	lai_temp = biomass_to_lai(biomass(ipts,ivm,ileaf,icarbon), ivm)
2657
2658	! This converts LAI to a height. I took this from functional_allocation.
2659	max_height = lai_temp * lai_to_height(ivm)
2660
2661	! The units here are per square meter, so the volume of the cube is
2662	! just the height11.
2663	IF (max_height .GT. min_stomate) THEN
2664	FAVD(ipts,ivm)=lai_temp / max_height
2665	ELSE
2666	FAVD(ipts,ivm)=zero
2667	ENDIF
2668	ENDIF
2669
2670
2671	IF (bavard.GE.2) WRITE(numout,*) 'Leaving calculate_favd'
2672
2673	END SUBROUTINE calculate_favd
2674
2675	!! ================================================================================================================================
2676	!! FUNCTION : calculate_laieff_fit
2677	!!
2678	!>\BRIEF :
2679	!!
2680	!! DESCRIPTION :
2681	!!
2682	!! RECENT CHANGE(S): None
2683	!!
2684	!! RETURN VALUE : calculate_laieff_fit
2685	!!
2686	!! REFERENCE(S) :
2687	!!
2688	!! FLOWCHART : None
2689	!! \n
2690	!_ ================================================================================================================================
2691
2692	REAL(r_std) FUNCTION calculate_laieff_fit(test_angle, laieff_fit)
2693	!! 0 Variable and parameter declaration
2694
2695	!! 0.1 Input variables
2696	REAL(r_std),INTENT(IN) :: test_angle !! the cosine of the solar angle
2697	TYPE(laieff_type),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
2698
2699	!! 0.2 Output variables
2700
2701	!! 0.3 Modified variables
2702
2703	!! 0.4 Local variables
2704	!_ ================================================================================================================================
2705
2706	! y = a + bx + cx*2 + dx*3 + ex**4
2707	calculate_laieff_fit=laieff_fit%a + test_angle * laieff_fit%b + &
2708	test_angle*2 laieff_fit%c + test_angle*3 laieff_fit%d + &
2709	test_angle*4 laieff_fit%e
2710
2711
2712
2713	END FUNCTION calculate_laieff_fit
2714
2715	END MODULE stomate_laieff

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

Context Navigation

source: branches/publications/ORCHIDEE_CAN_r3069/src_stomate/stomate_laieff.f90 @ 7475

Download in other formats: