Context Navigation

sapiens_product_use.f90 @ 7475

Last change on this file since 7475 was 2442, checked in by sebastiaan.luyssaert, 10 years ago
DEV: replaced most STOPs by ipslerr. This version compiles but has not been tested
File size: 55.2 KB

Line
1	! =================================================================================================================================
2	! MODULE : sapiens_product_use
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 Following lateral transport, calculate C-stored in the product pools
10	!! and the CO2-release from product decomposition
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! REFERENCE(S) : None
17	!!
18	!! SVN :
19	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/stomate_lcchange.f90 $
20	!! $Date: 2013-07-03 10:16:15 +0200 (Wed, 03 Jul 2013) $
21	!! $Revision: 1356 $
22	!! \n
23	!_ ================================================================================================================================
24
25	MODULE sapiens_product_use
26
27	! modules used:
28
29	USE stomate_data
30	USE pft_parameters
31	USE constantes
32
33	IMPLICIT NONE
34
35	PRIVATE
36	PUBLIC basic_product_use, basic_product_init, dim_product_use
37
38	CONTAINS
39
40
41	!! ================================================================================================================================
42	!! SUBROUTINE : basic_product_use
43	!!
44	!>\BRIEF Calculate the C stored in biomass-based products and the
45	!! CO2 release through their decomposition
46	!!
47	!! DESCRIPTION : The basic approach distinguished 3 pools (as introduced by
48	!! Shilong Piao based on Houghton. This module does NOT make use of the
49	!! dimensions of the woody harvest or of the management and cutting flags.
50	!! All wood is distributed across the short, medium and long product pools.
51	!! This basically means that wood use in ..., 1600, 1700, 1800, 1900 and 2000
52	!! was identical.\n
53	!!
54	!! RECENT CHANGE(S) : None
55	!!
56	!! MAIN OUTPUT VARIABLE(S) : prod_s, prod_m, prod_l, flux_s, flux_m, flux_l,
57	!! cflux_prod_s cflux_prod_m and cflux_prod1_l
58	!!
59	!! REFERENCES : None
60	!!
61	!! FLOWCHART : None
62	!!
63	!_ ================================================================================================================================
64
65	SUBROUTINE basic_product_use(npts, dt_days, harvest_pool_bound, harvest_pool, &
66	harvest_type, harvest_cut, harvest_area, &
67	cflux_prod_s, cflux_prod_m, cflux_prod_l, prod_s, &
68	prod_m, prod_l, prod_s_total, prod_m_total, &
69	prod_l_total, cflux_prod_total, flux_s, flux_m, &
70	flux_l, resolution)
71
72	IMPLICIT NONE
73
74
75	!! 0. Variable and parameter declaration
76
77	!! 0.1 Input variables
78
79	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
80	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
81	!! (days)
82	REAL(r_std),DIMENSION(:,:), INTENT(in) :: resolution !! The length in m of each grid-box in X and Y
83	REAL(r_std), DIMENSION(:), INTENT(in) :: harvest_pool_bound!! The boundaries of the diameter classes
84	!! in the wood harvest pools @tex $(m)$ @endtex
85
86	!! 0.2 Output variables
87
88	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_s !! C-released during first years (short term)
89	!! following land cover change
90	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
91	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_m !! Total annual release from decomposition of
92	!! the medium-lived product pool
93	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
94	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_l !! Total annual release from decomposition of
95	!! the long-lived product pool
96	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
97	REAL(r_std), DIMENSION(:), INTENT(out) :: prod_s_total !! Total products remaining in the short-lived
98	!! pool after the annual decomposition
99	!! @tex $(gC m^{-2})$ @endtex
100	REAL(r_std), DIMENSION(:), INTENT(out) :: prod_m_total !! Total products remaining in the medium-lived
101	!! pool after the annual decomposition
102	!! @tex $(gC m^{-2})$ @endtex
103	REAL(r_std), DIMENSION(:), INTENT(out) :: prod_l_total !! Total products remaining in the long-lived
104	!! pool after the annual release
105	!! @tex $(gC m^{-2})$ @endtex
106	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_total !! Total flux from decomposition of the short,
107	!! medium and long lived product pools
108	!! @tex $(gC m^{-1} dt^{-1})$ @endtex
109
110	!! 0.3 Modified variables
111
112	REAL(r_std), DIMENSION(:,0:), INTENT(inout) :: prod_s !! Short-lived product pool after the annual
113	!! release of each compartment (short + 1 :
114	!! input from year of land cover change)
115	!! @tex ($gC pixel^{-1}$) @endtex
116	REAL(r_std), DIMENSION(:,0:), INTENT(inout) :: prod_m !! Medium-lived product pool after the annual
117	!! release of each compartment (medium + 1 :
118	!! input from year of land cover change)
119	!! @tex ($gC pixel^{-1}$) @endtex
120	REAL(r_std), DIMENSION(:,0:), INTENT(inout) :: prod_l !! long-lived product pool after the annual
121	!! release of each compartment (long + 1 :
122	!! input from year of land cover change)
123	!! @tex ($gC pixel^{-1}$) @endtex
124	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_s !! Annual release from the short-lived product
125	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
126	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_m !! Annual release from the medium-lived product
127	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
128	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_l !! Annual release from the long-lived product
129	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
130	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: harvest_pool !! The wood and biomass that have been
131	!! havested by humans @tex $(gC)$ @endtex
132	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_type !! Type of management that resulted
133	!! in the harvest (unitless)
134	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_cut !! Type of cutting that was used for the harvest
135	!! (unitless)
136	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_area !! Harvested area (m^{2})
137
138
139	!! 0.4 Local variables
140
141	INTEGER(i_std) :: ipts, ivm, iage !! Indices (unitless)
142	INTEGER(i_std) :: imbc !! Indices (unitless)
143	REAL(r_std) :: exported_biomass !! Total biomass exported from the PFT (gC)
144	REAL(r_std) :: residual !! Residual @tex $(gC m^{-1})$ @endtex
145	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
146	!! mass balance chech for this routine
147	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
148	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
149	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
150	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start pool of this routine
151	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
152	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! End pool of this routine
153	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
154	REAL(r_std), DIMENSION(npts) :: pixel_area !! surface area of the pixel @tex ($m^{2}$) @endtex
155	!_ ================================================================================================================================
156
157	IF (bavard.GE.3) WRITE(numout,*) 'Entering basic product use'
158
159
160	!! 1. initialization
161
162	!! 1.1 Fluxes and pools
163	prod_s(:,0) = zero
164	prod_m(:,0) = zero
165	prod_l(:,0) = zero
166	prod_s_total(:) = zero
167	prod_m_total(:) = zero
168	prod_l_total(:) = zero
169	cflux_prod_s(:) = zero
170	cflux_prod_m(:) = zero
171	cflux_prod_l(:) = zero
172	cflux_prod_total(:) = zero
173
174	!! 1.2 Pixel area
175	! Surface area (m2) of the pixel
176	pixel_area(:) = resolution(:,1) * resolution(:,2)
177
178	!! 1.3 Mass balance closure
179	! The biomass harvest pool shouldn't be multiplied by veget_max
180	! its units are already in gC pixel-1. The total amount for the
181	! pixel was stored in harvest_pool. Account for the harvest and
182	! wood product pools. There are no longer PFTs
183	pool_start(:,:,:) = zero
184	pool_start(:,1,icarbon) = &
185	( SUM(SUM(harvest_pool(:,:,:,icarbon),3),2) + &
186	SUM(prod_l(:,:),2) + SUM(prod_m(:,:),2) + &
187	SUM(prod_s(:,:),2) )/pixel_area(:)
188
189	!! 1.4 Logical check of coeff_lcchange
190	! Within a PFT the coeff_lcchange should add up to one. Else the
191	! carbon balance will not be closed! Do not include PFT 1 in this
192	! check because its coeff_lcchange were set to undef (-9999).
193	DO ivm = 2,nvm
194
195	IF ( (coeff_lcchange_s(ivm) + coeff_lcchange_m(ivm) + &
196	coeff_lcchange_l(ivm) .LT. un - min_stomate) .OR. &
197	(coeff_lcchange_s(ivm) + coeff_lcchange_m(ivm) + &
198	coeff_lcchange_l(ivm) .GT. un + min_stomate) ) THEN
199
200	WRITE(numout,*) 'Error: coeff_lcchange in sapiens_product_use.f90 ' //&
201	'do not add up to 1 in PFT - cannot close the mass balance, ',ivm
202	CALL ipslerr_p (3,'sapiens_product_use', &
203	'basic_product_use','coeff_lcchange do not add up to 1 ','')
204	ENDIF
205
206	ENDDO
207
208
209	!! 2. Partition harvest in a short, medium and long-lived product pool
210
211	DO ipts = 1,npts
212
213	!! 2.1 Calculate the initial mass of the product pool (gC)
214	DO ivm = 2,nvm
215
216	!+++NOTE+++
217	! For the moment all harvest is lumped together and we don't
218	! make use of harvest_type and harvest_cut to refine the product
219	! use. If these variables are to be used be aware that for the
220	! moment these values can be overwritten within a simulation.
221	! Although the sequence of the calls in stomate_lpj is believed
222	! to protect against this, it is not enforced by the code. For
223	! example, if we first thin and then harvest, wood of both
224	! actions will be stored in the harvest pool but it will be
225	! recored as just harvest (thinning will be overwritten by harvest)
226	exported_biomass = SUM(harvest_pool(ipts,ivm,:,icarbon))
227	!+++++++++++
228
229	! Attribute the harvested biomass to pools with a different
230	! turnover time. The turnover times are PFT depend which
231	! accounts for the different uses different PFT's are put to i.e.
232	! wood from forests vs grass or crops. The values stored in
233	! harvest_pool are absolute numbers gC (for that pixel), hence,
234	! the spatial extent of the LCC change, harvest or thinning has
235	! already been accounted for.
236	prod_s(ipts,0) = prod_s(ipts,0) + &
237	(coeff_lcchange_s(ivm) * exported_biomass)
238	prod_m(ipts,0) = prod_m(ipts,0) + &
239	(coeff_lcchange_m(ivm) * exported_biomass)
240	prod_l(ipts,0) = prod_l(ipts,0) + &
241	(coeff_lcchange_l(ivm) * exported_biomass)
242
243	ENDDO
244
245	!! 2.2 Check whether all the harvest was used as is to be expected
246	! The precision we are aiming for is 10-8 per m-2. The harvest pool
247	! is expressed for the whole pixel
248	residual = SUM(SUM(harvest_pool(ipts,:,:,icarbon),2)) - &
249	prod_s(ipts,0) - prod_m(ipts,0) - prod_l(ipts,0)
250
251	IF (residual/pixel_area(ipts) .GT. min_stomate .OR. &
252	residual/pixel_area(ipts) .LT. - min_stomate) THEN
253
254	WRITE (numout,*) 'Error: some harvest in sapiens_product_use has '//&
255	'not been attributed'
256	WRITE (numout,*) 'Residual in pixel, ', ipts, residual
257	CALL ipslerr_p (3,'sapiens_product_use', &
258	'basic_product_use','some harvest has not been attributed','')
259
260	ELSE
261
262	! We could still have residual within the precision at the
263	! m2 level but outside the precision at the pixel level.
264	! Move the residual to the short lived pool
265	prod_s(ipts,0) = prod_s(ipts,0) + residual
266
267	! All the carbon in the harvest pool was allocated to the product
268	! pools the harvest pool can be set to zero so we can us it with
269	! confidence in the mass balance calculation
270	harvest_pool(ipts,:,:,icarbon) = zero
271	harvest_type(ipts,:) = zero
272	harvest_cut(ipts,:) = zero
273	harvest_area(ipts,:) = zero
274
275	ENDIF
276
277
278	!! 2.3 Update the long-turnover pool content
279	! Update the long-turnover pool content following flux emission
280	! of a linear decay (1/longivety) of the initial carbon input.
281	! This must be implemented with a counter going from ::nlong to
282	! zero because this way it takes ::nlong years before the intial
283	! pool has been entirely consumed. If a more intuitive approach
284	! would have been used in which the counter goes from zero to
285	! ::nlong, the pool would be depleted in the first year.
286	! Update the long-turnover pool content following flux emission.
287	DO iage = nlong,2,-1
288
289	! flux and prod are large numbers. At the last time step 10-5
290	! gC remains which is outside the precision. Therefore, this
291	! carbon is moved to the flux
292	IF (iage .EQ. nlong) THEN
293
294	cflux_prod_l(ipts) = cflux_prod_l(ipts) + prod_l(ipts,iage-1)
295	flux_l(ipts,iage) = prod_l(ipts,iage-1)
296	prod_l(ipts,iage) = zero
297
298	ELSE
299
300	! ::flux_l has INTENT(inout), it thus has the value from the
301	! previuous year at the start of these calculations. Copy
302	! the flux from the previous year to this year and update
303	! the product pool
304	cflux_prod_l(ipts) = cflux_prod_l(ipts) + flux_l(ipts,iage)
305	prod_l(ipts,iage) = prod_l(ipts,iage-1) - flux_l(ipts,iage-1)
306	flux_l(ipts,iage) = flux_l(ipts,iage-1)
307
308	ENDIF
309
310	ENDDO
311
312	! Treat the first year separately because it is not included in
313	! the loop above
314	cflux_prod_l(ipts) = cflux_prod_l(ipts) + flux_l(ipts,1)
315
316	! Each year 1/::nlong of the initial pool is decomposed. This fixed
317	! amount of decomposition is moved to a higher age class each year.
318	! The units of ::flux_m are gC (absolute number as is the product
319	! pool). At this point in the calculation, the flux is NOT a flux
320	! yet but is also a pool
321	flux_l(ipts,1) = (un/nlong) * prod_l(ipts,0)
322	prod_l(ipts,1) = prod_l(ipts,0) - flux_l(ipts,1)
323	prod_l(ipts,0) = zero
324
325	!! 2.4 Update the medium-turnover pool content
326	! Update the medium-turnover pool content following flux emission.
327	! For comments see the section above. The concept is identical.
328	DO iage = nmedium,2,-1
329
330	IF (iage .EQ. nmedium) THEN
331
332	cflux_prod_m(ipts) = cflux_prod_m(ipts) + prod_m(ipts,iage-1)
333	flux_m(ipts,iage) = prod_m(ipts,iage-1)
334	prod_m(ipts,iage) = zero
335
336	ELSE
337
338	cflux_prod_m(ipts) = cflux_prod_m(ipts) + flux_m(ipts,iage)
339	prod_m(ipts,iage) = prod_m(ipts,iage-1) - flux_m(ipts,iage-1)
340	flux_m(ipts,iage) = flux_m(ipts,iage-1)
341
342	ENDIF
343
344	ENDDO
345
346	cflux_prod_m(ipts) = cflux_prod_m(ipts) + flux_m(ipts,1)
347	flux_m(ipts,1) = (un/nmedium) * prod_m(ipts,0)
348	prod_m(ipts,1) = prod_m(ipts,0) - flux_m(ipts,1)
349	prod_m(ipts,0) = zero
350
351	!! 2.5 Update the short-turnover pool
352	! The length of the short-turnover pool will most likely be
353	! less than length of the regular loop. The regular loop is
354	! now conditional.
355	! NOTE: this code has only been tested for nshort = 1 and
356	! nshort = 2
357	IF (nshort .GE. 4) THEN
358
359	DO iage = nshort,2,-1
360
361	IF (iage .EQ. nshort) THEN
362
363	cflux_prod_s(ipts) = cflux_prod_s(ipts) + prod_s(ipts,iage-1)
364	flux_s(ipts,iage) = prod_s(ipts,iage-1)
365	prod_s(ipts,iage) = zero
366
367	ELSE
368
369	cflux_prod_s(ipts) = cflux_prod_s(ipts) + flux_s(ipts,iage)
370	prod_s(ipts,iage) = prod_s(ipts,iage-1) - flux_s(ipts,iage-1)
371	flux_s(ipts,iage) = flux_s(ipts,iage-1)
372
373	ENDIF
374
375	ENDDO
376
377	cflux_prod_s(ipts) = cflux_prod_s(ipts) + flux_s(ipts,1)
378	flux_s(ipts,1) = (un/nshort) * prod_s(ipts,0)
379	prod_s(ipts,1) = prod_s(ipts,0) - flux_s(ipts,1)
380	prod_s(ipts,0) = zero
381
382	ELSEIF (nshort .EQ. 3) THEN
383
384	! Decompose the product pool from the last year
385	cflux_prod_s(ipts) = cflux_prod_s(ipts) + prod_s(ipts,2)
386	flux_s(ipts,3) = prod_s(ipts,2)
387	prod_s(ipts,3) = zero
388	cflux_prod_s(ipts) = cflux_prod_s(ipts) + flux_s(ipts,2)
389	prod_s(ipts,2) = prod_s(ipts,1) - flux_s(ipts,1)
390	flux_s(ipts,2) = flux_s(ipts,1)
391
392	! Calculate the linear product decomposition for this
393	! year's harvest
394	cflux_prod_s(ipts) = cflux_prod_s(ipts) + flux_s(ipts,1)
395	flux_s(ipts,1) = (un/nshort) * prod_s(ipts,0)
396	prod_s(ipts,1) = prod_s(ipts,0) - flux_s(ipts,1)
397	prod_s(ipts,0) = zero
398
399	ELSEIF (nshort .EQ. 2) THEN
400
401	! Decompose the product pool from last year
402	cflux_prod_s(ipts) = cflux_prod_s(ipts) + prod_s(ipts,1)
403	flux_s(ipts,2) = prod_s(ipts,1)
404	prod_s(ipts,2) = zero
405
406	! Calculate the linear product decomposition for this
407	! year's harvest
408	cflux_prod_s(ipts) = cflux_prod_s(ipts) + flux_s(ipts,1)
409	flux_s(ipts,1) = (un/nshort) * prod_s(ipts,0)
410	prod_s(ipts,1) = prod_s(ipts,0) - flux_s(ipts,1)
411	prod_s(ipts,0) = zero
412
413	ELSEIF (nshort .EQ. 1) THEN
414
415	! Everything decomposes in a single year. The short product
416	! pool remains empty.
417	cflux_prod_s(ipts) = cflux_prod_s(ipts) + prod_s(ipts,0)
418	flux_s(ipts,1) = prod_s(ipts,0)
419	prod_s(ipts,1) = zero
420	prod_s(ipts,0) = zero
421
422	ELSE
423
424	WRITE(numout,*) 'Error: flaw in the logic of IF-construct '//&
425	'(basic_product_use in sapiens_product_use.f90)'
426
427	ENDIF ! number of ages in short-turnover product pool
428
429	ENDDO ! #pixels - spatial domain
430
431
432	!! 3. Check mass balance closure
433
434	!! 3.1 Calculate pools at the end of the routine
435	! The biomass harvest pool shouldn't be multiplied by veget_max
436	! its in gC pixel-1 so divide by pixel_area to get a value in
437	! gC m-2. Account for the harvest and wood product pools. The harvest
438	! pool should be empty, it was added for completeness.
439	closure_intern = zero
440	pool_end = zero
441	pool_end(:,1,icarbon) = &
442	( SUM(SUM(harvest_pool(:,:,:,icarbon),3),2) + &
443	SUM(prod_l(:,:),2) + SUM(prod_m(:,:),2) + &
444	SUM(prod_s(:,:),2) ) / pixel_area(:)
445
446	!! 3.2 Calculate components of the mass balance
447	check_intern(:,1,iatm2land,icarbon) = zero
448	check_intern(:,1,iland2atm,icarbon) = -un * ( SUM(flux_s(:,:),2) + &
449	SUM(flux_m(:,:),2) + SUM(flux_l(:,:),2) ) / pixel_area(:)
450	check_intern(:,1,ilat2out,icarbon) = zero
451	check_intern(:,1,ilat2in,icarbon) = -un * zero
452	check_intern(:,1,ipoolchange,icarbon) = -un * (pool_end(:,1,icarbon) - &
453	pool_start(:,1,icarbon))
454	DO imbc = 1,nmbcomp
455	closure_intern(:,1,icarbon) = closure_intern(:,1,icarbon) + &
456	check_intern(:,1,imbc,icarbon)
457	ENDDO
458
459	! 3.3 Write outcome
460	DO ipts=1,npts
461	DO ivm=1,nvm
462	IF(ABS(closure_intern(ipts,ivm,icarbon)) .LE. min_stomate)THEN
463	IF (ld_massbal) WRITE(numout,*) 'Mass balance closure in basic_product_use'
464	ELSE
465	WRITE(numout,*) 'Error: mass balance is not closed in basic_product_use'
466	WRITE(numout,*) ' ipts,ivm; ', ipts,ivm
467	WRITE(numout,*) ' Difference is, ', closure_intern(ipts,ivm,icarbon)
468	WRITE(numout,*) ' pool_end,pool_start: ', pool_end(ipts,ivm,icarbon), pool_start(ipts,ivm,icarbon)
469	WRITE(numout,*) ' fluxes: ', flux_s(ipts,ivm),flux_m(ipts,ivm),flux_l(ipts,ivm)
470	WRITE(numout,*) ' pixel_area: ', pixel_area(ipts)
471	IF(ld_stop)THEN
472	CALL ipslerr_p (3,'basic_product_use', 'Mass balance error.','','')
473	ENDIF
474	ENDIF
475	ENDDO
476	ENDDO
477
478	!! 4. Convert the pools into fluxes and aggregate
479
480	! Convert pools into fluxes
481	! The pools are given in the absolute amount of carbon for the whole pixel
482	! and the harvest was accumulated over the year. The unit is thus
483	! gc pixel-1 year-1. These fluxes should be expressed as gC m-2 dt-1. Therefore,
484	! divide by the resolution of the pixel and the number of time steps within
485	! a year
486	cflux_prod_s(:) = cflux_prod_s(:) / one_year * dt_days / pixel_area(:)
487	cflux_prod_m(:) = cflux_prod_m(:) / one_year * dt_days / pixel_area(:)
488	cflux_prod_l(:) = cflux_prod_l(:) / one_year * dt_days / pixel_area(:)
489
490	! Calculate aggregate values
491	! cflux_prod_total is now in gC m-2 dt-1. The prod_x_total are in gC m-2
492	cflux_prod_total(:) = cflux_prod_s(:) + cflux_prod_m(:) + cflux_prod_l(:)
493	prod_s_total(:) = SUM(prod_s(:,:),2) / pixel_area(:)
494	prod_m_total(:) = SUM(prod_m(:,:),2) / pixel_area(:)
495	prod_l_total(:) = SUM(prod_l(:,:),2) / pixel_area(:)
496
497
498	IF (bavard.GE.4) WRITE(numout,*) 'Leaving basic product use'
499
500
501	END SUBROUTINE basic_product_use
502
503
504	!! ================================================================================================================================
505	!! SUBROUTINE : dim_product_use
506	!!
507	!>\BRIEF Calculate the C stored in biomass-based products and the
508	!! CO2 release through their decomposition.
509	!!
510	!! DESCRIPTION : The basic approach distinguished 3 pools (as introduced by
511	!! Shilong Piao based on Houghton. Wood harvest with small dimensions is being
512	!! used as fire wood the remaining wood goes into the short, medium and long
513	!! product pools. This implies that when the dimensions of the harvest change,
514	!! the product pools will change as well.\n
515	!!
516	!! RECENT CHANGE(S) : None
517	!!
518	!! MAIN OUTPUT VARIABLE(S) : prod_s, prod_m, prod_l, flux_s, flux_m, flux_l,
519	!! cflux_prod_s cflux_prod_m and cflux_prod1_l
520	!!
521	!! REFERENCES : None
522	!!
523	!! FLOWCHART : None
524	!!
525	!_ ================================================================================================================================
526
527	SUBROUTINE dim_product_use(npts, dt_days, harvest_pool_bound, harvest_pool, &
528	harvest_type, harvest_cut, harvest_area, &
529	cflux_prod_s, cflux_prod_m, cflux_prod_l, prod_s, &
530	prod_m, prod_l, prod_s_total, prod_m_total, &
531	prod_l_total, cflux_prod_total, flux_s, flux_m, &
532	flux_l, resolution)
533
534	IMPLICIT NONE
535
536
537	!! 0. Variable and parameter declaration
538
539	!! 0.1 Input variables
540
541	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
542	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
543	!! (days)
544	REAL(r_std),DIMENSION(:,:), INTENT(in) :: resolution !! The length in m of each grid-box in X and Y
545	REAL(r_std), DIMENSION(:), INTENT(in) :: harvest_pool_bound!! The boundaries of the diameter classes
546	!! in the wood harvest pools @tex $(m)$ @endtex
547
548	!! 0.2 Output variables
549
550	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_s !! C-released during first years (short term)
551	!! following land cover change
552	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
553	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_m !! Total annual release from decomposition of
554	!! the medium-lived product pool
555	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
556	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_l !! Total annual release from decomposition of
557	!! the long-lived product pool
558	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
559	REAL(r_std), DIMENSION(:), INTENT(out) :: prod_s_total !! Total products remaining in the short-lived
560	!! pool after the annual decomposition
561	!! @tex $(gC m^{-2})$ @endtex
562	REAL(r_std), DIMENSION(:), INTENT(out) :: prod_m_total !! Total products remaining in the medium-lived
563	!! pool after the annual decomposition
564	!! @tex $(gC m^{-2})$ @endtex
565	REAL(r_std), DIMENSION(:), INTENT(out) :: prod_l_total !! Total products remaining in the long-lived
566	!! pool after the annual release
567	!! @tex $(gC m^{-2})$ @endtex
568	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_total !! Total flux from decomposition of the short,
569	!! medium and long lived product pools
570	!! @tex $(gC m^{-1} dt^{-1})$ @endtex
571
572	!! 0.3 Modified variables
573
574	REAL(r_std), DIMENSION(:,0:), INTENT(inout) :: prod_s !! Short-lived product pool after the annual
575	!! release of each compartment (short + 1 :
576	!! input from year of land cover change)
577	!! @tex ($gC pixel^{-1}$) @endtex
578	REAL(r_std), DIMENSION(:,0:), INTENT(inout) :: prod_m !! Medium-lived product pool after the annual
579	!! release of each compartment (medium + 1 :
580	!! input from year of land cover change)
581	!! @tex ($gC pixel^{-1}$) @endtex
582	REAL(r_std), DIMENSION(:,0:), INTENT(inout) :: prod_l !! long-lived product pool after the annual
583	!! release of each compartment (long + 1 :
584	!! input from year of land cover change)
585	!! @tex ($gC pixel^{-1}$) @endtex
586	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_s !! Annual release from the short-lived product
587	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
588	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_m !! Annual release from the medium-lived product
589	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
590	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_l !! Annual release from the long-lived product
591	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
592	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: harvest_pool !! The wood and biomass that have been
593	!! havested by humans @tex $(gC)$ @endtex
594	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_type !! Type of management that resulted
595	!! in the harvest (unitless)
596	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_cut !! Type of cutting that was used for the harvest
597	!! (unitless)
598	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_area !! Harvested area (m^{2})
599
600
601	!! 0.4 Local variables
602
603	INTEGER(i_std) :: ipts, ivm, iage !! Indices (unitless)
604	INTEGER(i_std) :: imbc, idia !! Indices (unitless)
605	REAL(r_std) :: exported_biomass !! Total biomass exported from the PFT (gC)
606	REAL(r_std) :: fuelwood_biomass !! Harvest with small dimensions that will
607	!! be used as fire wood (gC)
608	REAL(r_std) :: timber_biomass !! Harvest with larger dimensions that will
609	!! be used for all kinds of uses (gC)
610	REAL(r_std) :: residual !! Residual @tex $(gC m^{-1})$ @endtex
611	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
612	!! mass balance chech for this routine
613	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
614	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
615	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
616	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start pool of this routine
617	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
618	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! End pool of this routine
619	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
620	REAL(r_std), DIMENSION(npts) :: pixel_area !! surface area of the pixel @tex ($m^{2}$) @endtex
621	!_ ================================================================================================================================
622
623	IF (bavard.GE.3) WRITE(numout,*) 'Entering dimensional product use'
624
625
626	!! 1. initialization
627
628	!! 1.1 Fluxes and pools
629	prod_s(:,0) = zero
630	prod_m(:,0) = zero
631	prod_l(:,0) = zero
632	prod_s_total(:) = zero
633	prod_m_total(:) = zero
634	prod_l_total(:) = zero
635	cflux_prod_s(:) = zero
636	cflux_prod_m(:) = zero
637	cflux_prod_l(:) = zero
638	cflux_prod_total(:) = zero
639
640	!! 1.2 Pixel area
641	! Surface area (m2) of the pixel
642	pixel_area(:) = resolution(:,1) * resolution(:,2)
643
644	!! 1.3 Mass balance closure
645	! The biomass harvest pool shouldn't be multiplied by veget_max
646	! its units are already in gC pixel-1. The total amount for the
647	! pixel was stored in harvest_pool. Account for the harvest and
648	! wood product pools. There are no longer PFTs
649	pool_start(:,:,:) = zero
650	pool_start(:,1,icarbon) = &
651	( SUM(SUM(harvest_pool(:,:,:,icarbon),3),2) + &
652	SUM(prod_l(:,:),2) + SUM(prod_m(:,:),2) + &
653	SUM(prod_s(:,:),2) )/pixel_area(:)
654
655	!! 1.4 Logical check of coeff_lcchange
656	! Within a PFT the coeff_lcchange should add up to one. Else the
657	! carbon balance will not be closed! Do not include PFT 1 in this
658	! check because its coeff_lcchange were set to undef (-9999).
659	DO ivm = 2,nvm
660
661	IF ( (coeff_lcchange_s(ivm) + coeff_lcchange_m(ivm) + &
662	coeff_lcchange_l(ivm) .LT. un - min_stomate) .OR. &
663	(coeff_lcchange_s(ivm) + coeff_lcchange_m(ivm) + &
664	coeff_lcchange_l(ivm) .GT. un + min_stomate) ) THEN
665
666	WRITE(numout,*) 'Error: coeff_lcchange in sapiens_product_use.f90 ' //&
667	'do not add up to 1 in PFT - cannot close the mass balance, ',ivm
668	CALL ipslerr_p (3,'sapiens_product_use', &
669	'dim_product_use','coeff_lcchange do not add up to 1 ','')
670
671	ENDIF
672
673	ENDDO
674
675
676	!! 2. Partition harvest in a short, medium and long-lived product pool
677
678	DO ipts = 1,npts
679
680	!! 2.1 Calculate the initial mass of the product pool (gC)
681	DO ivm = 2,nvm
682
683	!+++NOTE+++
684	! For the moment all harvest is lumped together and we don't
685	! make use of harvest_type and harvest_cut to refine the product
686	! use. If these variables are to be used be aware that for the
687	! moment these values can be overwritten within a simulation.
688	! Although the sequence of the calls in stomate_lpj is believed
689	! to protect against this, it is not enforced by the code. For
690	! example, if we first thin and then harvest, wood of both
691	! actions will be stored in the harvest pool but it will be
692	! recored as just harvest (thinning will be overwritten by
693	! harvest and no warning will be given)
694	!++++++++++
695	fuelwood_biomass = zero
696	timber_biomass = zero
697
698	! Grasses and crops
699	IF (fuelwood_diameter(ivm) .LT. min_stomate) THEN
700
701	!Grasses and crops do not provide fuelwood. They do have
702	! a very short turnover time. The values stored in
703	! harvest_pool are absolute numbers gC (for that pixel), hence,
704	! the spatial extent of the LCC change, harvest or thinning has
705	! already been accounted for.
706	exported_biomass = SUM(harvest_pool(ipts,ivm,:,icarbon))
707	prod_s(ipts,0) = prod_s(ipts,0) + &
708	(coeff_lcchange_s(ivm) * exported_biomass)
709	prod_m(ipts,0) = prod_m(ipts,0) + &
710	(coeff_lcchange_m(ivm) * exported_biomass)
711	prod_l(ipts,0) = prod_l(ipts,0) + &
712	(coeff_lcchange_l(ivm) * exported_biomass)
713
714	! Forests
715	ELSE
716
717	! This is a woody PFT, check the diameter of the harvested wood
718	timber_biomass = SUM(harvest_pool(ipts,ivm,:,icarbon))
719	DO idia = 1,ndia_harvest
720
721	IF (harvest_pool_bound(idia) .LT. fuelwood_diameter(ivm)) THEN
722
723	! If the diameter of the harvest pool is below the treshold
724	! the wood is used as fuelwood and will end up in the short
725	! lived pool
726	fuelwood_biomass = fuelwood_biomass + &
727	harvest_pool(ipts,ivm,idia,icarbon)
728
729	ENDIF
730
731	ENDDO ! harvest diameter
732
733	! Calculate the long lived pools as the residual of all
734	! wood mines the fuel wood.
735	timber_biomass = timber_biomass - fuelwood_biomass
736
737	! Attribute the harvested biomass to pools with a different
738	! turnover time. The turnover times are PFT depend which
739	! accounts for the different uses different PFT's are put to i.e.
740	! wood from forests vs grass or crops. The values stored in
741	! harvest_pool are absolute numbers gC (for that pixel), hence,
742	! the spatial extent of the LCC change, harvest or thinning has
743	! already been accounted for.
744	prod_s(ipts,0) = prod_s(ipts,0) + fuelwood_biomass
745	IF (coeff_lcchange_m(ivm) + coeff_lcchange_l(ivm) .GT. &
746	min_stomate) THEN
747
748	! If the condition is satisfied, calculate the fractions
749	! if it is not satisfied there is no medium and long lived
750	! pool for the products. Nothing should then be done. Avoid
751	! divide by zero
752	prod_m(ipts,0) = prod_m(ipts,0) + &
753	(coeff_lcchange_m(ivm)/ &
754	(coeff_lcchange_m(ivm)+coeff_lcchange_l(ivm)) * &
755	timber_biomass)
756	prod_l(ipts,0) = prod_l(ipts,0) + &
757	(coeff_lcchange_l(ivm)/ &
758	(coeff_lcchange_m(ivm)+coeff_lcchange_l(ivm)) * &
759	timber_biomass)
760
761	ENDIF ! Check coeff_lcchange
762
763	ENDIF ! grasses/crops vs forests
764
765	ENDDO ! loop over PFTs
766
767	!! 2.2 Check whether all the harvest was used as is to be expected
768	! The precision we are aiming for is 10-8 per m-2. The harvest pool
769	! is expressed for the whole pixel
770	residual = SUM(SUM(harvest_pool(ipts,:,:,icarbon),2)) - &
771	prod_s(ipts,0) - prod_m(ipts,0) - prod_l(ipts,0)
772
773	IF (residual/pixel_area(ipts) .GT. min_stomate .OR. &
774	residual/pixel_area(ipts) .LT. - min_stomate) THEN
775
776	WRITE (numout,*) 'Error: some harvest in sapiens_product_use has '//&
777	'not been attributed'
778	WRITE (numout,*) 'Residual in pixel, ', ipts, residual
779	CALL ipslerr_p (3,'sapiens_product_use', &
780	'basic_product_use','coeff_lcchange do not add up to 1 ','')
781
782	ELSE
783
784	! We could still have residual within the precision at the
785	! m2 level but outside the precision at the pixel level.
786	! Move the residual to the short lived pool
787	prod_s(ipts,0) = prod_s(ipts,0) + residual
788
789	! All the carbon in the harvest pool was allocated to the product
790	! pools the harvest pool can be set to zero so we can us it with
791	! confidence in the mass balance calculation
792	harvest_pool(ipts,:,:,icarbon) = zero
793	harvest_type(ipts,:) = zero
794	harvest_cut(ipts,:) = zero
795	harvest_area(ipts,:) = zero
796
797	ENDIF
798
799
800	!! 2.3 Update the long-turnover pool content
801	! Update the long-turnover pool content following flux emission
802	! of a linear decay (1/longivety) of the initial carbon input.
803	! This must be implemented with a counter going from ::nlong to
804	! zero because this way it takes ::nlong years before the intial
805	! pool has been entirely consumed. If a more intuitive approach
806	! would have been used in which the counter goes from zero to
807	! ::nlong, the pool would be depleted in the first year.
808	! Update the long-turnover pool content following flux emission.
809	DO iage = nlong,2,-1
810
811	! flux and prod are large numbers. At the last time step 10-5
812	! gC remains which is outside the precision. Therefore, this
813	! carbon is moved to the flux
814	IF (iage .EQ. nlong) THEN
815
816	cflux_prod_l(ipts) = cflux_prod_l(ipts) + prod_l(ipts,iage-1)
817	flux_l(ipts,iage) = prod_l(ipts,iage-1)
818	prod_l(ipts,iage) = zero
819
820	ELSE
821
822	! ::flux_l has INTENT(inout), it thus has the value from the
823	! previuous year at the start of these calculations. Copy
824	! the flux from the previous year to this year and update
825	! the product pool
826	cflux_prod_l(ipts) = cflux_prod_l(ipts) + flux_l(ipts,iage)
827	prod_l(ipts,iage) = prod_l(ipts,iage-1) - flux_l(ipts,iage-1)
828	flux_l(ipts,iage) = flux_l(ipts,iage-1)
829
830	ENDIF
831
832	ENDDO
833
834	! Treat the first year separately because it is not included in
835	! the loop above
836	cflux_prod_l(ipts) = cflux_prod_l(ipts) + flux_l(ipts,1)
837
838	! Each year 1/::nlong of the initial pool is decomposed. This fixed
839	! amount of decomposition is moved to a higher age class each year.
840	! The units of ::flux_m are gC (absolute number as is the product
841	! pool). At this point in the calculation, the flux is NOT a flux
842	! yet but is also a pool
843	flux_l(ipts,1) = (un/nlong) * prod_l(ipts,0)
844	prod_l(ipts,1) = prod_l(ipts,0) - flux_l(ipts,1)
845	prod_l(ipts,0) = zero
846
847	!! 2.4 Update the medium-turnover pool content
848	! Update the medium-turnover pool content following flux emission.
849	! For comments see the section above. The concept is identical.
850	DO iage = nmedium,2,-1
851
852	IF (iage .EQ. nmedium) THEN
853
854	cflux_prod_m(ipts) = cflux_prod_m(ipts) + prod_m(ipts,iage-1)
855	flux_m(ipts,iage) = prod_m(ipts,iage-1)
856	prod_m(ipts,iage) = zero
857
858	ELSE
859
860	cflux_prod_m(ipts) = cflux_prod_m(ipts) + flux_m(ipts,iage)
861	prod_m(ipts,iage) = prod_m(ipts,iage-1) - flux_m(ipts,iage-1)
862	flux_m(ipts,iage) = flux_m(ipts,iage-1)
863
864	ENDIF
865
866	ENDDO
867
868	cflux_prod_m(ipts) = cflux_prod_m(ipts) + flux_m(ipts,1)
869	flux_m(ipts,1) = (un/nmedium) * prod_m(ipts,0)
870	prod_m(ipts,1) = prod_m(ipts,0) - flux_m(ipts,1)
871	prod_m(ipts,0) = zero
872
873	!! 2.5 Update the short-turnover pool
874	! The length of the short-turnover pool will most likely be
875	! less than length of the regular loop. The regular loop is
876	! now conditional.
877	! NOTE: this code has only been tested for nshort = 1 and
878	! nshort = 2
879	IF (nshort .GE. 4) THEN
880
881	DO iage = nshort,2,-1
882
883	IF (iage .EQ. nshort) THEN
884
885	cflux_prod_s(ipts) = cflux_prod_s(ipts) + prod_s(ipts,iage-1)
886	flux_s(ipts,iage) = prod_s(ipts,iage-1)
887	prod_s(ipts,iage) = zero
888
889	ELSE
890
891	cflux_prod_s(ipts) = cflux_prod_s(ipts) + flux_s(ipts,iage)
892	prod_s(ipts,iage) = prod_s(ipts,iage-1) - flux_s(ipts,iage-1)
893	flux_s(ipts,iage) = flux_s(ipts,iage-1)
894
895	ENDIF
896
897	ENDDO
898
899	cflux_prod_s(ipts) = cflux_prod_s(ipts) + flux_s(ipts,1)
900	flux_s(ipts,1) = (un/nshort) * prod_s(ipts,0)
901	prod_s(ipts,1) = prod_s(ipts,0) - flux_s(ipts,1)
902	prod_s(ipts,0) = zero
903
904	ELSEIF (nshort .EQ. 3) THEN
905
906	! Decompose the product pool from the last year
907	cflux_prod_s(ipts) = cflux_prod_s(ipts) + prod_s(ipts,2)
908	flux_s(ipts,3) = prod_s(ipts,2)
909	prod_s(ipts,3) = zero
910	cflux_prod_s(ipts) = cflux_prod_s(ipts) + flux_s(ipts,2)
911	prod_s(ipts,2) = prod_s(ipts,1) - flux_s(ipts,1)
912	flux_s(ipts,2) = flux_s(ipts,1)
913
914	! Calculate the linear product decomposition for this
915	! year's harvest
916	cflux_prod_s(ipts) = cflux_prod_s(ipts) + flux_s(ipts,1)
917	flux_s(ipts,1) = (un/nshort) * prod_s(ipts,0)
918	prod_s(ipts,1) = prod_s(ipts,0) - flux_s(ipts,1)
919	prod_s(ipts,0) = zero
920
921	ELSEIF (nshort .EQ. 2) THEN
922
923	! Decompose the product pool from last year
924	cflux_prod_s(ipts) = cflux_prod_s(ipts) + prod_s(ipts,1)
925	flux_s(ipts,2) = prod_s(ipts,1)
926	prod_s(ipts,2) = zero
927
928	! Calculate the linear product decomposition for this
929	! year's harvest
930	cflux_prod_s(ipts) = cflux_prod_s(ipts) + flux_s(ipts,1)
931	flux_s(ipts,1) = (un/nshort) * prod_s(ipts,0)
932	prod_s(ipts,1) = prod_s(ipts,0) - flux_s(ipts,1)
933	prod_s(ipts,0) = zero
934
935	ELSEIF (nshort .EQ. 1) THEN
936
937	! Everything decomposes in a single year. The short product
938	! pool remains empty.
939	cflux_prod_s(ipts) = cflux_prod_s(ipts) + prod_s(ipts,0)
940	flux_s(ipts,1) = prod_s(ipts,0)
941	prod_s(ipts,1) = zero
942	prod_s(ipts,0) = zero
943
944	ELSE
945
946	WRITE(numout,*) 'Error: flaw in the logic of IF-construct '//&
947	'(basic_product_use in sapiens_product_use.f90)'
948
949
950	ENDIF ! number of ages in short-turnover product pool
951
952	ENDDO ! #pixels - spatial domain
953
954
955	!! 3. Check mass balance closure
956
957	!! 3.1 Calculate pools at the end of the routine
958	! The biomass harvest pool shouldn't be multiplied by veget_max
959	! its in gC pixel-1 so divide by pixel_area to get a value in
960	! gC m-2. Account for the harvest and wood product pools. The harvest
961	! pool should be empty, it was added for completeness.
962	closure_intern = zero
963	pool_end = zero
964	pool_end(:,1,icarbon) = &
965	( SUM(SUM(harvest_pool(:,:,:,icarbon),3),2) + &
966	SUM(prod_l(:,:),2) + SUM(prod_m(:,:),2) + &
967	SUM(prod_s(:,:),2) ) / pixel_area(:)
968
969	!! 3.2 Calculate components of the mass balance
970	check_intern(:,1,iatm2land,icarbon) = zero
971	check_intern(:,1,iland2atm,icarbon) = -un * ( SUM(flux_s(:,:),2) + &
972	SUM(flux_m(:,:),2) + SUM(flux_l(:,:),2) ) / pixel_area(:)
973	check_intern(:,1,ilat2out,icarbon) = zero
974	check_intern(:,1,ilat2in,icarbon) = -un * zero
975	check_intern(:,1,ipoolchange,icarbon) = -un * (pool_end(:,1,icarbon) - &
976	pool_start(:,1,icarbon))
977	DO imbc = 1,nmbcomp
978	closure_intern(:,1,icarbon) = closure_intern(:,1,icarbon) + &
979	check_intern(:,1,imbc,icarbon)
980	ENDDO
981
982	! 3.3 Write outcome
983	DO ipts=1,npts
984	DO ivm=1,nvm
985	IF(ABS(closure_intern(ipts,ivm,icarbon)) .LE. min_stomate)THEN
986	IF (ld_massbal) WRITE(numout,*) 'Mass balance closure in basic_product_use'
987	ELSE
988	WRITE(numout,*) 'Error: mass balance is not closed in basic_product_use'
989	WRITE(numout,*) ' ipts,ivm; ', ipts,ivm
990	WRITE(numout,*) ' Difference is, ', closure_intern(ipts,ivm,icarbon)
991	WRITE(numout,*) ' pool_end,pool_start: ', pool_end(ipts,ivm,icarbon), pool_start(ipts,ivm,icarbon)
992	WRITE(numout,*) ' fluxes: ', flux_s(ipts,ivm),flux_m(ipts,ivm),flux_l(ipts,ivm)
993	WRITE(numout,*) ' pixel_area: ', pixel_area(ipts)
994	IF(ld_stop)THEN
995	CALL ipslerr_p (3,'dim_product_use', 'Mass balance error.','','')
996	ENDIF
997	ENDIF
998	ENDDO
999	ENDDO
1000
1001	!! 4. Convert the pools into fluxes and aggregate
1002
1003	! Convert pools into fluxes
1004	! The pools are given in the absolute amount of carbon for the whole pixel
1005	! and the harvest was accumulated over the year. The unit is thus
1006	! gc pixel-1 year-1. These fluxes should be expressed as gC m-2 dt-1. Therefore,
1007	! divide by the resolution of the pixel and the number of time steps within
1008	! a year
1009	cflux_prod_s(:) = cflux_prod_s(:) / one_year * dt_days / pixel_area(:)
1010	cflux_prod_m(:) = cflux_prod_m(:) / one_year * dt_days / pixel_area(:)
1011	cflux_prod_l(:) = cflux_prod_l(:) / one_year * dt_days / pixel_area(:)
1012
1013	! Calculate aggregate values
1014	! cflux_prod_total is now in gC m-2 dt-1. The prod_x_total are in gC m-2
1015	cflux_prod_total(:) = cflux_prod_s(:) + cflux_prod_m(:) + cflux_prod_l(:)
1016	prod_s_total(:) = SUM(prod_s(:,:),2) / pixel_area(:)
1017	prod_m_total(:) = SUM(prod_m(:,:),2) / pixel_area(:)
1018	prod_l_total(:) = SUM(prod_l(:,:),2) / pixel_area(:)
1019
1020
1021	IF (bavard.GE.4) WRITE(numout,*) 'Leaving basic product use'
1022
1023
1024	END SUBROUTINE dim_product_use
1025
1026	!! ================================================================================================================================
1027	!! SUBROUTINE : basic_product_init
1028	!!
1029	!>\BRIEF Initialize the product use variables when product use is not
1030	!! calculated
1031	!!
1032	!! DESCRIPTION : Initialize the product use variables when product use is not
1033	!! calculated
1034	!!
1035	!! RECENT CHANGE(S) : None
1036	!!
1037	!! MAIN OUTPUT VARIABLE(S) : prod_s, prod_m, prod_l, flux_s, flux_m, flux_l,
1038	!! cflux_prod_s cflux_prod_m and cflux_prod1_l
1039	!!
1040	!! REFERENCES : None
1041	!!
1042	!! FLOWCHART : None
1043	!!
1044	!_ ================================================================================================================================
1045
1046	SUBROUTINE basic_product_init(cflux_prod_s, cflux_prod_m, cflux_prod_l, &
1047	prod_s_total, prod_m_total, prod_l_total, flux_s, &
1048	flux_m, flux_l, cflux_prod_total)
1049
1050	IMPLICIT NONE
1051
1052
1053	!! 0. Variable and parameter declaration
1054
1055	!! 0.1 Input variables
1056
1057	!! 0.2 Output variables
1058
1059	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_s !! C-released during first years (short term)
1060	!! following land cover change
1061	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
1062	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_m !! Total annual release from decomposition of
1063	!! the medium-lived product pool
1064	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
1065	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_l !! Total annual release from decomposition of
1066	!! the long-lived product pool
1067	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
1068	REAL(r_std), DIMENSION(:), INTENT(out) :: prod_s_total !! Total products remaining in the short-lived
1069	!! pool after the annual decomposition
1070	!! @tex $(gC m^{-2})$ @endtex
1071	REAL(r_std), DIMENSION(:), INTENT(out) :: prod_m_total !! Total products remaining in the medium-lived
1072	!! pool after the annual decomposition
1073	!! @tex $(gC m^{-2})$ @endtex
1074	REAL(r_std), DIMENSION(:), INTENT(out) :: prod_l_total !! Total products remaining in the long-lived
1075	!! pool after the annual release
1076	!! @tex $(gC m^{-2})$ @endtex
1077	REAL(r_std), DIMENSION(:), INTENT(out) :: cflux_prod_total !! Total flux from decomposition of the short,
1078	!! medium and long lived product pools
1079	!! @tex $(gC m^{-1} dt^{-1})$ @endtex
1080
1081	!! 0.3 Modified variables
1082
1083	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_s !! Annual release from the short-lived product
1084	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
1085	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_m !! Annual release from the medium-lived product
1086	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
1087	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_l !! Annual release from the long-lived product
1088	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
1089
1090	!! 0.4 Local variables
1091
1092	!_ ================================================================================================================================
1093
1094	IF (bavard.GE.3) WRITE(numout,*) 'Entering basic product init'
1095
1096	!! 1. initialization
1097
1098	cflux_prod_s(:) = zero
1099	cflux_prod_m(:) = zero
1100	cflux_prod_l(:) = zero
1101	prod_s_total(:) = zero
1102	prod_m_total(:) = zero
1103	prod_l_total(:) = zero
1104	flux_s(:,:) = zero
1105	flux_m(:,:) = zero
1106	flux_l(:,:) = zero
1107	cflux_prod_total(:) = zero
1108
1109	IF (bavard.GE.4) WRITE(numout,*) 'Leaving basic product init'
1110
1111	END SUBROUTINE basic_product_init
1112
1113	END MODULE sapiens_product_use

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

Context Navigation

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

Download in other formats: