Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_stomate/sapiens_product_use.f90 @ 7852

Last change on this file since 7852 was 7317, checked in by sebastiaan.luyssaert, 3 years ago
Bug fixes for product pools. Contributes to ticket #774. Changes proposed by Chao
File size: 65.4 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	USE grid
33	USE function_library, ONLY: check_mass_balance
34
35	IMPLICIT NONE
36
37	PRIVATE
38	PUBLIC basic_product_use, product_init, dim_product_use
39
40	INTEGER(i_std), SAVE :: printlev_loc !! Local level of text output for this module
41	!$OMP THREADPRIVATE(printlev_loc)
42	LOGICAL, SAVE :: firstcall_sapiens_product_use = .TRUE. !! first call flag
43	!$OMP THREADPRIVATE(firstcall_sapiens_product_use)
44
45	CONTAINS
46
47
48	!! ================================================================================================================================
49	!! SUBROUTINE : basic_product_use
50	!!
51	!>\BRIEF Calculate the C stored in biomass-based products and the
52	!! CO2 release through their decomposition
53	!!
54	!! DESCRIPTION : The basic approach distinguished 3 pools (as introduced by
55	!! Shilong Piao based on Houghton. This module does NOT make use of the
56	!! dimensions of the woody harvest or of the management and cutting flags.
57	!! All wood is distributed across the short, medium and long product pools.
58	!! This basically means that wood use in ..., 1600, 1700, 1800, 1900 and 2000
59	!! was identical. Wood product pools of medium and long turnover times start
60	!! decomposition at the next year of harvest; while short-residence wood
61	!! product pool starts decomposition at the year of harvest.\n
62	!!
63	!! RECENT CHANGE(S) : None
64	!!
65	!! MAIN OUTPUT VARIABLE(S) : prod_s, prod_m, prod_l, flux_s, flux_m, flux_l,
66	!! flux_prod_s flux_prod_m and flux_prod_l
67	!!
68	!! REFERENCES : None
69	!!
70	!! FLOWCHART : None
71	!!
72	!_ ================================================================================================================================
73
74	SUBROUTINE basic_product_use(npts, dt_days, harvest_pool_bound, harvest_pool_acc, &
75	harvest_type, harvest_cut, harvest_area, &
76	flux_prod_s, flux_prod_m, flux_prod_l, prod_s, &
77	prod_m, prod_l, prod_s_total, prod_m_total, &
78	prod_l_total, flux_prod_total, flux_s, flux_m, &
79	flux_l, veget_max, flux_s_pft)
80
81	IMPLICIT NONE
82
83
84	!! 0. Variable and parameter declaration
85
86	!! 0.1 Input variables
87
88	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
89	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
90	!! (days)
91	REAL(r_std), DIMENSION(:), INTENT(in) :: harvest_pool_bound!! The boundaries of the diameter classes
92	!! in the wood harvest pools @tex $(m)$ @endtex
93	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! Passed so we can check_mass_balance
94	!! Making it OPTIONAL would have been nicer.
95
96
97	!! 0.2 Output variables
98
99	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_s !! C-released during first years (short term)
100	!! following land cover change
101	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
102	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_m !! Total annual release from decomposition of
103	!! the medium-lived product pool
104	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
105	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_l !! Total annual release from decomposition of
106	!! the long-lived product pool
107	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
108	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: prod_s_total !! Total products remaining in the short-lived
109	!! pool after the annual decomposition
110	!! @tex $(gC)$ @endtex
111	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: prod_m_total !! Total products remaining in the medium-lived
112	!! pool after the annual decomposition
113	!! @tex $(gC)$ @endtex
114	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: prod_l_total !! Total products remaining in the long-lived
115	!! pool after the annual release
116	!! @tex $(gC)$ @endtex
117	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_total !! Total flux from decomposition of the short,
118	!! medium and long lived product pools
119	!! @tex $(gC m^{-1} dt^{-1})$ @endtex
120	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_s_pft !! Total flux from decomposition the same year of the lcc and harvest
121	!! @tex $(gC m^{-1} dt^{-1})$ @endtex
122
123
124	!! 0.3 Modified variables
125
126	REAL(r_std), DIMENSION(:,0:,:,:,:), INTENT(inout) :: prod_s !! Short-lived product pool after the annual
127	!! release of each compartment (short + 1 :
128	!! input from year of land cover change)
129	!! @tex ($gC pixel^{-1}$) @endtex
130	REAL(r_std), DIMENSION(:,0:,:,:,:), INTENT(inout) :: prod_m !! Medium-lived product pool after the annual
131	!! release of each compartment (medium + 1 :
132	!! input from year of land cover change)
133	!! @tex ($gC pixel^{-1}$) @endtex
134	REAL(r_std), DIMENSION(:,0:,:,:,:), INTENT(inout) :: prod_l !! long-lived product pool after the annual
135	!! release of each compartment (long + 1 :
136	!! input from year of land cover change)
137	!! @tex ($gC pixel^{-1}$) @endtex
138	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: flux_s !! Annual release from the short-lived product
139	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
140	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: flux_m !! Annual release from the medium-lived product
141	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
142	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: flux_l !! Annual release from the long-lived product
143	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
144	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: harvest_pool_acc !! The wood and biomass that have been
145	!! havested by humans @tex $(gC pixel-1)$ @endtex
146	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_type !! Type of management that resulted
147	!! in the harvest (unitless)
148	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_cut !! Type of cutting that was used for the harvest
149	!! (unitless)
150	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_area !! Harvested area (m^{2})
151
152
153	!! 0.4 Local variables
154
155	INTEGER(i_std) :: ipts, ivm, iage !! Indices (unitless)
156	INTEGER(i_std) :: imbc, iele, ilan !! Indices (unitless)
157	INTEGER(i_std) :: ilct !! Indices (unitless)
158	REAL(r_std), DIMENSION(nelements,nlctypes) :: harvest_acc_tmp !! Accumulated harvest per land cover type (gC pixel-1)
159	REAL(r_std), DIMENSION(nelements,nlanduse) :: exported_biomass !! Exported biomass per land use type from the PFT (gC pixel-1)
160	REAL(r_std), DIMENSION(nelements) :: residual !! Residual @tex $(gC m^{-1})$ @endtex
161	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
162	!! mass balance chech for this routine
163	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
164	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
165	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
166	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start pool of this routine
167	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
168	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! End pool of this routine
169	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
170
171	!_ ================================================================================================================================
172
173
174	!! 1. initialization
175
176	!! 1.1 Initialize local printlev only first time one subroutine in the module is called
177	IF (firstcall_sapiens_product_use) THEN
178	printlev_loc=get_printlev('sapiens_product_use')
179	firstcall_sapiens_product_use=.FALSE.
180	END IF
181
182	IF (printlev_loc.GE.2) WRITE(numout,*) 'Entering basic product use'
183
184	!! 1.2 Fluxes and pools
185	prod_s(:,0,:,:,:) = zero
186	prod_m(:,0,:,:,:) = zero
187	prod_l(:,0,:,:,:) = zero
188	prod_s_total(:,:,:,:) = zero
189	prod_m_total(:,:,:,:) = zero
190	prod_l_total(:,:,:,:) = zero
191	flux_prod_s(:,:,:,:) = zero
192	flux_prod_m(:,:,:,:) = zero
193	flux_prod_l(:,:,:,:) = zero
194	flux_prod_total(:,:,:,:) = zero
195	flux_s_pft(:,:,:,:) = zero
196
197	!! 1.4 Initialize check for mass balance closure
198	IF (err_act.GT.1) THEN
199
200	! The biomass harvest pool shouldn't be multiplied by veget_max
201	! its units are already in gC pixel-1. The total amount for the
202	! pixel was stored in harvest_pool_acc. Account for the harvest and
203	! wood product pools. There are no longer PFTs
204	check_intern(:,:,:,:) = zero
205	pool_start(:,:,:) = zero
206	DO iele = 1,nelements
207	pool_start(:,1,iele) = &
208	( SUM(SUM(harvest_pool_acc(:,:,:,iele),3),2) + &
209	SUM(SUM(SUM(prod_l(:,:,iele,:,:),2),2),2) + &
210	SUM(SUM(SUM(prod_m(:,:,iele,:,:),2),2),2) + &
211	SUM(SUM(SUM(prod_s(:,:,iele,:,:),2),2),2) ) / area(:)
212	ENDDO
213
214	ENDIF ! err_act.GT.1
215
216	!! 1.5 Logical check of coeff_lcchange
217	! Within a PFT the coeff_lcchange should add up to one. Else the
218	! carbon balance will not be closed! Do not include PFT 1 in this
219	! check because its coeff_lcchange were set to undef (-9999).
220	DO ivm = 2,nvm
221
222	IF ( ABS(coeff_lcchange_s(ivm) + coeff_lcchange_m(ivm) + &
223	coeff_lcchange_l(ivm)-un) .GT. min_stomate ) THEN
224
225	WRITE(numout,*) 'Error: coeff_lcchange in sapiens_product_use.f90 ' //&
226	'do not add up to 1 in PFT - cannot close the mass balance, ',ivm
227	CALL ipslerr_p (3,'sapiens_product_use', &
228	'basic_product_use','coeff_lcchange do not add up to 1 ','')
229	ENDIF
230
231	ENDDO
232
233
234	!! 2. Partition harvest in a short, medium and long-lived product pool
235
236	DO ipts = 1,npts
237
238	harvest_acc_tmp(:,:) = zero
239
240	!! 2.1 Calculate the initial mass of the product pool (gC)
241	DO ivm = 2,nvm
242
243	! Determine the land cover type
244	IF (is_tree(ivm)) THEN
245	ilct = iforest
246	ELSEIF (.NOT. is_tree(ivm) .AND. natural(ivm)) THEN
247	ilct = igrass
248	ELSEIF (.NOT. is_tree(ivm) .AND. .NOT. natural(ivm)) THEN
249	ilct = icrop
250	ELSE
251	CALL ipslerr_p (3,'sapiens_product_use.f90',&
252	'not clear to which land cover type the PFT belongs',&
253	'check the code', '')
254	END IF
255
256	harvest_acc_tmp(:,ilct) = harvest_acc_tmp(:,ilct) + &
257	SUM(harvest_pool_acc(ipts,ivm,:,:),1)
258
259	! Product use is calculated separatly for the biomass export due
260	! to lcc and harvest. Be aware that for the
261	! moment these values can be overwritten within a simulation.
262	! Although the sequence of the calls in stomate_lpj is believed
263	! to protect against this, it is not enforced by the code. For
264	! example, if we first thin and then harvest, wood of both
265	! actions will be stored in the harvest pool but it will be
266	! recored as just harvest (thinning will be overwritten by harvest)
267	exported_biomass(:,:) = zero
268
269	IF(harvest_cut(ipts,ivm).EQ.icut_lcc_wood .OR. &
270	harvest_cut(ipts,ivm).EQ.icut_lcc_res) THEN
271
272	! All biomass exported following lcc.
273	exported_biomass(:,ilcc) = SUM(harvest_pool_acc(ipts,ivm,:,:),1)
274
275
276	ELSE
277
278	! All biomass exported following some sort of harvest
279	exported_biomass(:,iharvest) = SUM(harvest_pool_acc(ipts,ivm,:,:),1)
280
281	ENDIF
282
283	!! 2.2 Process the exported biomass
284	! Attribute the harvested biomass to pools with a different
285	! turnover time. The turnover times are PFT depend which
286	! accounts for the different uses different PFT's are put to i.e.
287	! wood from forests vs grass or crops. The values stored in
288	! harvest_pool_acc are absolute numbers gC (for that pixel), hence,
289	! the spatial extent of the LCC change, harvest or thinning has
290	! already been accounted for.
291	prod_s(ipts,0,:,:,ilct) = prod_s(ipts,0,:,:,ilct) + &
292	(coeff_lcchange_s(ivm) * exported_biomass(:,:))
293	prod_m(ipts,0,:,:,ilct) = prod_m(ipts,0,:,:,ilct) + &
294	(coeff_lcchange_m(ivm) * exported_biomass(:,:))
295	prod_l(ipts,0,:,:,ilct) = prod_l(ipts,0,:,:,ilct) + &
296	(coeff_lcchange_l(ivm) * exported_biomass(:,:))
297
298	! Calculate how much C and N is released back into the
299	! atmosphere the same year (needed for AR6 output)
300	! Note: when nshort=1, this makes completely sense. For nshort>1,
301	! this only accounts for the carbon release happening within the
302	! the year of harvest.
303	flux_s_pft(ipts,ivm,:,ilct) = SUM(harvest_pool_acc(ipts,ivm,:,:),dim=1) * &
304	coeff_lcchange_s(ivm) / nshort
305
306	END DO ! loop over PFTs
307
308	DO ilct = 1,nlctypes
309
310	! Check whether all the harvest was used as is to be expected
311	! The precision we are aiming for is 10-8 per m-2. The harvest pool
312	! is expressed for the whole pixel
313	residual(:) = harvest_acc_tmp(:,ilct) - &
314	SUM(prod_s(ipts,0,:,:,ilct),2) - &
315	SUM(prod_m(ipts,0,:,:,ilct),2) - &
316	SUM(prod_l(ipts,0,:,:,ilct),2)
317
318	DO iele = 1,nelements
319
320	IF (ABS(residual(iele)/area(ipts)) .GT. min_stomate) THEN
321
322	WRITE (numout,*) 'Error: some harvest in '//&
323	'sapiens_product_use has not been attributed'
324	WRITE (numout,*) 'Residual in pixel, iele, ', &
325	iele, ipts, residual
326	CALL ipslerr_p (3,'sapiens_product_use', &
327	'basic_product_use','some harvest has not been attributed',&
328	'')
329
330	ELSE
331
332	! We could still have residuals within the precision at the
333	! m2 level but outside the precision at the pixel level.
334	! Move the residual to the short lived pool.
335	IF ( prod_s(ipts,0,iele,iharvest,ilct).GT.zero) THEN
336
337	! If there was harvest put the residual in the product pool of
338	! wood harvest
339	prod_s(ipts,0,iele,iharvest,ilct) = prod_s(ipts,0,iele,iharvest,ilct) + &
340	residual(iele)
341
342	ELSE
343
344	! There was no harvest but there is some residual so the residual
345	! should come from land cover change
346	prod_s(ipts,0,iele,ilcc,ilct) = prod_s(ipts,0,iele,ilcc,ilct) + &
347	residual(iele)
348
349	END IF
350
351	END IF
352
353	END DO ! iele
354
355	END DO ! land cover types
356
357	! All the carbon in the harvest pool was allocated to the product
358	! pools the harvest pool can be set to zero so we can us it with
359	! confidence in the mass balance calculation
360	harvest_pool_acc(ipts,:,:,:) = zero
361	harvest_type(ipts,:) = zero
362	harvest_cut(ipts,:) = zero
363	harvest_area(ipts,:) = zero
364
365	!! 2.3 Update the long-turnover pool content
366	! Update the long-turnover pool content following flux emission
367	! of a linear decay (1/longivety) of the initial carbon input.
368	! This must be implemented with a counter going from ::nlong to
369	! zero because this way it takes ::nlong years before the intial
370	! pool has been entirely consumed. If a more intuitive approach
371	! would have been used in which the counter goes from zero to
372	! ::nlong, the pool would be depleted in the first year.
373	! Update the long-turnover pool content following flux emission.
374	DO iage = nlong,2,-1
375
376	! below codes are a simple bookkeeping one, note that flux_l contains
377	! the constant annual value that needs to be removed every year.
378	! The key is therefore to just shift their positions
379
380	! Suppose nlong = 10; logic of the following lines are:
381	!
382	! flux_prod = flux_prod + flux(10)
383	! prod(10) = prod(9) - flux(9) !shift the 9th to 10th
384	! flux(10) = flux(9) ! this is to prepare for next year
385
386	! flux_prod = flux_prod + flux(9)
387	! prod(9) = prod(8) - flux(8) !shift the 8th to 9th
388	! flux(9) = flux(8)
389
390	! ...
391
392	! flux_prod = flux_prod + flux(2)
393	! prod(2) = prod(1) - flux(1) !shift the 1st to 2nd
394	! flux(2) = flux(1)
395
396	flux_prod_l(ipts,:,:,:) = flux_prod_l(ipts,:,:,:) + flux_l(ipts,iage,:,:,:)
397	prod_l(ipts,iage,:,:,:) = prod_l(ipts,iage-1,:,:,:) - flux_l(ipts,iage-1,:,:,:)
398	flux_l(ipts,iage,:,:,:) = flux_l(ipts,iage-1,:,:,:)
399
400	ENDDO
401
402	! Treat the first year separately
403	flux_prod_l(ipts,:,:,:) = flux_prod_l(ipts,:,:,:) + flux_l(ipts,1,:,:,:)
404
405	! Each year 1/::nlong of the initial pool is decomposed. This fixed
406	! amount of decomposition is moved to a higher age class each year.
407	! The units of ::flux_m are gC (absolute number as is the product
408	! pool). At this point in the calculation, the flux is NOT a flux
409	! yet but is also a pool. Note harvest pool generated this year only
410	! starts to decompose the next year.
411	flux_l(ipts,1,:,:,:) = (un/nlong) * prod_l(ipts,0,:,:,:)
412	prod_l(ipts,1,:,:,:) = prod_l(ipts,0,:,:,:)
413	prod_l(ipts,0,:,:,:) = zero
414
415	!! 2.4 Update the medium-turnover pool content
416	! Update the medium-turnover pool content following flux emission.
417	! For comments see the section above. The concept is identical.
418	DO iage = nmedium,2,-1
419	flux_prod_m(ipts,:,:,:) = flux_prod_m(ipts,:,:,:) + flux_m(ipts,iage,:,:,:)
420	prod_m(ipts,iage,:,:,:) = prod_m(ipts,iage-1,:,:,:) - flux_m(ipts,iage-1,:,:,:)
421	flux_m(ipts,iage,:,:,:) = flux_m(ipts,iage-1,:,:,:)
422	ENDDO
423
424	flux_prod_m(ipts,:,:,:) = flux_prod_m(ipts,:,:,:) + flux_m(ipts,1,:,:,:)
425	flux_m(ipts,1,:,:,:) = (un/nmedium) * prod_m(ipts,0,:,:,:)
426	prod_m(ipts,1,:,:,:) = prod_m(ipts,0,:,:,:)
427	prod_m(ipts,0,:,:,:) = zero
428
429	!! 2.5 Update the short-turnover pool
430	! The length of the short-turnover pool will most likely be
431	! less than length of the regular loop. Different from long- and
432	! medium-residence pool, we suppose decomposition of short pool
433	! start from the year of harvest, rather than the next year of
434	! harvest year.
435	IF (nshort .GE. 3) THEN
436
437	! note we start from nshort-1 here,because starting from the harvest
438	! year, it takes 'nshort-1' years for the pool to completely return to
439	! the atmosphere.
440	DO iage = nshort-1,2,-1
441	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + flux_s(ipts,iage,:,:,:)
442	prod_s(ipts,iage,:,:,:) = prod_s(ipts,iage-1,:,:,:) - flux_s(ipts,iage-1,:,:,:)
443	flux_s(ipts,iage,:,:,:) = flux_s(ipts,iage-1,:,:,:)
444	ENDDO
445
446	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + flux_s(ipts,1,:,:,:)
447	flux_s(ipts,1,:,:,:) = (un/nshort) * prod_s(ipts,0,:,:,:)
448	! this line accounts for the decomposition happening at the year
449	! of harvest.
450	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + flux_s(ipts,1,:,:,:)
451	prod_s(ipts,1,:,:,:) = prod_s(ipts,0,:,:,:) - flux_s(ipts,1,:,:,:)
452	prod_s(ipts,0,:,:,:) = zero
453
454	ELSEIF (nshort .EQ. 2) THEN
455	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + flux_s(ipts,1,:,:,:)
456	flux_s(ipts,1,:,:,:) = (un/nshort) * prod_s(ipts,0,:,:,:)
457	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + flux_s(ipts,1,:,:,:)
458	prod_s(ipts,1,:,:,:) = prod_s(ipts,0,:,:,:) - flux_s(ipts,1,:,:,:)
459	prod_s(ipts,0,:,:,:) = zero
460
461	ELSEIF (nshort .EQ. 1) THEN
462
463	! Everything decomposes in a single year. The short product
464	! pool remains empty.
465	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + prod_s(ipts,0,:,:,:)
466	prod_s(ipts,0,:,:,:) = zero
467
468	ELSE
469
470	WRITE(numout,*) 'Error: flaw in the logic of IF-construct '//&
471	'(basic_product_use in sapiens_product_use.f90)'
472	CALL ipslerr_p(3,'ERROR: product_use',&
473	'Logic flaw in if-construct','','')
474
475	ENDIF ! number of ages in short-turnover product pool
476
477	ENDDO ! #pixels - spatial domain
478
479
480	!! 3. Check numerical consistency of this routine
481
482	IF (err_act.GT.1) THEN
483
484	! 3.1 Mass balance closure
485	! 3.1.1 Calculate final biomass
486	! The biomass harvest pool shouldn't be multiplied by veget_max
487	! its in gC pixel-1 so divide by area to get a value in
488	! gC m-2. Account for the harvest and wood product pools. The harvest
489	! pool should be empty, it was added for completeness.
490	pool_end = zero
491	DO iele=1,nelements
492	pool_end(:,1,iele) = &
493	( SUM(SUM(harvest_pool_acc(:,:,:,iele),3),2) + &
494	SUM(SUM(SUM(prod_l(:,:,iele,:,:),2),2),2) + &
495	SUM(SUM(SUM(prod_m(:,:,iele,:,:),2),2),2) + &
496	SUM(SUM(SUM(prod_s(:,:,iele,:,:),2),2),2) ) / area(:)
497	ENDDO
498
499	!! 3.1.2 Calculate mass balance
500	! Common processes
501	DO iele=1,nelements
502	check_intern(:,1,iland2atm,iele) = &
503	-un * ( SUM(SUM(flux_prod_s(:,iele,:,:) + &
504	flux_prod_m(:,iele,:,:) + &
505	flux_prod_l(:,iele,:,:),2),2) ) / area(:)
506	check_intern(:,1,ipoolchange,iele) = &
507	-un * (pool_end(:,1,iele) - &
508	pool_start(:,1,iele))
509	ENDDO
510
511	closure_intern = zero
512	DO imbc = 1,nmbcomp
513	DO iele=1,nelements
514	closure_intern(:,1,iele) = closure_intern(:,1,iele) + &
515	check_intern(:,1,imbc,iele)
516	ENDDO
517	ENDDO
518
519	! 3.1.3 Check mass balance closure
520	CALL check_mass_balance("basic_product_use", closure_intern, npts, &
521	pool_end, pool_start, veget_max, 'products')
522
523	ENDIF ! err_act.GT.1
524
525	!! 4. Convert the pools into fluxes and aggregate
526
527	! Convert pools into fluxes
528	! The pools are given in the absolute amount of carbon for the
529	! whole pixel and the harvest was accumulated over the year.
530	! The unit is thus gc pixel-1 year-1. These fluxes should be
531	! expressed as gC m-2 dt-1. Therefore, divide by the area
532	! of the pixel and the number of time steps within a year
533	DO iele = 1,nelements
534
535	DO ilan = 1,nlanduse
536
537	DO ilct = 1,nlctypes
538
539	! Convert flux_prod_total and flux_prod_x to gC m-2 dt-1. Calculate
540	! flux_prod_l as the residual for better mass conservation.
541	flux_prod_total(:,iele,ilan,ilct) = (flux_prod_s(:,iele,ilan,ilct) + &
542	flux_prod_m(:,iele,ilan,ilct) + flux_prod_l(:,iele,ilan,ilct)) / area(:)
543	flux_prod_s(:,iele,ilan,ilct) = flux_prod_s(:,iele,ilan,ilct) / area(:)
544	flux_prod_m(:,iele,ilan,ilct) = flux_prod_m(:,iele,ilan,ilct) / area(:)
545	flux_prod_l(:,iele,ilan,ilct) = flux_prod_total(:,iele,ilan,ilct) - &
546	flux_prod_s(:,iele,ilan,ilct) - flux_prod_m(:,iele,ilan,ilct)
547
548	! Convert prod_x_total in gC m-2. prod_x_total is not used in any mass
549	! balance check so there is no need to use a residual term to achieve
550	! higher precision.
551	prod_s_total(:,iele,ilan,ilct) = SUM(prod_s(:,:,iele,ilan,ilct),2)
552	prod_m_total(:,iele,ilan,ilct) = SUM(prod_m(:,:,iele,ilan,ilct),2)
553	prod_l_total(:,iele,ilan,ilct) = SUM(prod_l(:,:,iele,ilan,ilct),2)
554
555	END DO
556
557	END DO
558
559	! AR6 output
560	flux_s_pft(:,:,iele,:) = flux_s_pft(:,:,iele,:) / one_year * dt_days
561
562	END DO
563
564
565	IF(printlev.GE.4) WRITE(numout,*) 'Leaving basic product use'
566
567	END SUBROUTINE basic_product_use
568
569	!! ================================================================================================================================
570	!! SUBROUTINE : dim_product_use
571	!!
572	!>\BRIEF Calculate the C stored in biomass-based products and the
573	!! CO2 release through their decomposition.
574	!!
575	!! DESCRIPTION : The basic approach distinguished 3 pools (as introduced by
576	!! Shilong Piao based on Houghton. Wood harvest with small dimensions is being
577	!! used as fire wood the remaining wood goes into the short, medium and long
578	!! product pools. This implies that when the dimensions of the harvest change,
579	!! the product pools will change as well.\n
580	!!
581	!! RECENT CHANGE(S) : None
582	!!
583	!! MAIN OUTPUT VARIABLE(S) : prod_s, prod_m, prod_l, flux_s, flux_m, flux_l,
584	!! flux_prod_s flux_prod_m and flux_prod_l
585	!!
586	!! REFERENCES : None
587	!!
588	!! FLOWCHART : None
589	!!
590	!_ ================================================================================================================================
591
592	SUBROUTINE dim_product_use(npts, dt_days, harvest_pool_bound, harvest_pool_acc, &
593	harvest_type, harvest_cut, harvest_area, &
594	flux_prod_s, flux_prod_m, flux_prod_l, prod_s, &
595	prod_m, prod_l, prod_s_total, prod_m_total, &
596	prod_l_total, flux_prod_total, flux_s, flux_m, &
597	flux_l, veget_max, flux_s_pft)
598
599	IMPLICIT NONE
600
601
602	!! 0. Variable and parameter declaration
603
604	!! 0.1 Input variables
605
606	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
607	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
608	!! (days)
609	REAL(r_std), DIMENSION(:), INTENT(in) :: harvest_pool_bound!! The boundaries of the diameter classes
610	!! in the wood harvest pools @tex $(m)$ @endtex
611	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! Passed so we can check_mass_balance
612	!! Making it OPTIONAL would have been nicer.
613
614
615	!! 0.2 Output variables
616
617	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_s !! C-released during first years (short term)
618	!! following land cover change
619	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
620	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_m !! Total annual release from decomposition of
621	!! the medium-lived product pool
622	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
623	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_l !! Total annual release from decomposition of
624	!! the long-lived product pool
625	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
626	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: prod_s_total !! Total products remaining in the short-lived
627	!! pool after the annual decomposition
628	!! @tex $(gC)$ @endtex
629	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: prod_m_total !! Total products remaining in the medium-lived
630	!! pool after the annual decomposition
631	!! @tex $(gC)$ @endtex
632	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: prod_l_total !! Total products remaining in the long-lived
633	!! pool after the annual release
634	!! @tex $(gC)$ @endtex
635	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_total !! Total flux from decomposition of the short,
636	!! medium and long lived product pools
637	!! @tex $(gC m^{-1} dt^{-1})$ @endtex
638	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_s_pft !! Total flux from decomposition the same year of the lcc and harvest
639	!! @tex $(gC m^{-1} dt^{-1})$ @endtex
640
641	!! 0.3 Modified variables
642
643	REAL(r_std), DIMENSION(:,0:,:,:,:), INTENT(inout) :: prod_s !! Short-lived product pool after the annual
644	!! release of each compartment (short + 1 :
645	!! input from year of land cover change)
646	!! @tex ($gC pixel^{-1}$) @endtex
647	REAL(r_std), DIMENSION(:,0:,:,:,:), INTENT(inout) :: prod_m !! Medium-lived product pool after the annual
648	!! release of each compartment (medium + 1 :
649	!! input from year of land cover change)
650	!! @tex ($gC pixel^{-1}$) @endtex
651	REAL(r_std), DIMENSION(:,0:,:,:,:), INTENT(inout) :: prod_l !! long-lived product pool after the annual
652	!! release of each compartment (long + 1 :
653	!! input from year of land cover change)
654	!! @tex ($gC pixel^{-1}$) @endtex
655	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: flux_s !! Annual release from the short-lived product
656	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
657	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: flux_m !! Annual release from the medium-lived product
658	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
659	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: flux_l !! Annual release from the long-lived product
660	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
661	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: harvest_pool_acc !! The wood and biomass that have been
662	!! havested by humans @tex $(gC)$ @endtex
663	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_type !! Type of management that resulted
664	!! in the harvest (unitless)
665	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_cut !! Type of cutting that was used for the harvest
666	!! (unitless)
667	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_area !! Harvested area (m^{2})
668
669
670	!! 0.4 Local variables
671
672	INTEGER(i_std) :: ipts, ivm, iage !! Indices (unitless)
673	INTEGER(i_std) :: imbc, idia, iele !! Indices (unitless)
674	INTEGER(i_std) :: ilan, ilct !! Indices (unitless)
675	REAL(r_std), DIMENSION(nelements,nlctypes) :: harvest_acc_tmp !! Accumulated harvest per land cover type (gC pixel-1)
676	REAL(r_std), DIMENSION(nelements,nlanduse) :: exported_biomass !! Total biomass exported from the PFT (gC)
677	REAL(r_std), DIMENSION(nelements,nlanduse) :: fuelwood_biomass !! Harvest with small dimensions that will
678	!! be used as fire wood (gC)
679	REAL(r_std), DIMENSION(nelements,nlanduse) :: timber_biomass !! Harvest with larger dimensions that will
680	!! be used for all kinds of uses (gC)
681	REAL(r_std), DIMENSION(nelements) :: residual !! Residual @tex $(gC m^{-1})$ @endtex
682	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
683	!! mass balance chech for this routine
684	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
685	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
686	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
687	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start pool of this routine
688	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
689	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! End pool of this routine
690	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
691
692	!_ ================================================================================================================================
693
694
695	!! 1. initialization
696
697	!! 1.1 Initialize local printlev only first time one subroutine in the module is called
698	IF (firstcall_sapiens_product_use) THEN
699	printlev_loc=get_printlev('sapiens_product_use')
700	firstcall_sapiens_product_use=.FALSE.
701	END IF
702
703	IF (printlev_loc.GE.2) WRITE(numout,*) 'Entering dimensional product use'
704
705	!! 1.2 Fluxes and pools
706	prod_s(:,0,:,:,:) = zero
707	prod_m(:,0,:,:,:) = zero
708	prod_l(:,0,:,:,:) = zero
709	prod_s_total(:,:,:,:) = zero
710	prod_m_total(:,:,:,:) = zero
711	prod_l_total(:,:,:,:) = zero
712	flux_prod_s(:,:,:,:) = zero
713	flux_prod_m(:,:,:,:) = zero
714	flux_prod_l(:,:,:,:) = zero
715	flux_prod_total(:,:,:,:) = zero
716	flux_s_pft(:,:,:,:) = zero
717
718	!! 1.4 Initialize check for mass balance closure
719	IF (err_act.GT.1) THEN
720
721	! The biomass harvest pool shouldn't be multiplied by veget_max
722	! its units are already in gC pixel-1. The total amount for the
723	! pixel was stored in harvest_pool_acc. Account for the harvest and
724	! wood product pools. There are no longer PFTs
725	check_intern(:,:,:,:) = zero
726	pool_start(:,:,:) = zero
727	DO iele = 1,nelements
728	pool_start(:,1,iele) = &
729	( SUM(SUM(harvest_pool_acc(:,:,:,iele),3),2) + &
730	SUM(SUM(SUM(prod_l(:,:,iele,:,:),2),2),2) + &
731	SUM(SUM(SUM(prod_m(:,:,iele,:,:),2),2),2) + &
732	SUM(SUM(SUM(prod_s(:,:,iele,:,:),2),2),2) ) / area(:)
733	ENDDO
734
735	ENDIF ! err_act.GT.1
736
737	!! 1.4 Logical check of coeff_lcchange
738	! Within a PFT the coeff_lcchange should add up to one. Else the
739	! carbon balance will not be closed! Do not include PFT 1 in this
740	! check because its coeff_lcchange were set to undef (-9999).
741	DO ivm = 2,nvm
742
743	IF ( (coeff_lcchange_s(ivm) + coeff_lcchange_m(ivm) + &
744	coeff_lcchange_l(ivm) .LT. un - min_stomate) .OR. &
745	(coeff_lcchange_s(ivm) + coeff_lcchange_m(ivm) + &
746	coeff_lcchange_l(ivm) .GT. un + min_stomate) ) THEN
747
748	WRITE(numout,*) 'Error: coeff_lcchange in sapiens_product_use.f90 ' //&
749	'do not add up to 1 in PFT - cannot close the mass balance, ',ivm
750	CALL ipslerr_p (3,'sapiens_product_use', &
751	'dim_product_use','coeff_lcchange do not add up to 1 ','')
752
753	ENDIF
754
755	ENDDO
756
757
758	!! 2. Partition harvest in a short, medium and long-lived product pool
759
760	DO ipts = 1,npts
761
762	harvest_acc_tmp(:,:) = zero
763
764	!! 2.1 Calculate the initial mass of the product pool (gC)
765	DO ivm = 2,nvm
766
767	! Determine the land cover type
768	IF (is_tree(ivm)) THEN
769	ilct = iforest
770	ELSEIF (.NOT. is_tree(ivm) .AND. natural(ivm)) THEN
771	ilct = igrass
772	ELSEIF (.NOT. is_tree(ivm) .AND. .NOT. natural(ivm)) THEN
773	ilct = icrop
774	ELSE
775	CALL ipslerr_p (3,'sapiens_product_use.f90',&
776	'not clear to which land cover type the PFT belongs',&
777	'check the code', '')
778	END IF
779
780	harvest_acc_tmp(:,ilct) = harvest_acc_tmp(:,ilct) + &
781	SUM(harvest_pool_acc(ipts,ivm,:,:),1)
782
783	! Product use is calculated separatly for the biomass export due
784	! to lcc and harvest. Be aware that for the
785	! moment these values can be overwritten within a simulation.
786	! Although the sequence of the calls in stomate_lpj is believed
787	! to protect against this, it is not enforced by the code. For
788	! example, if we first thin and then harvest, wood of both
789	! actions will be stored in the harvest pool but it will be
790	! recored as just harvest (thinning will be overwritten by harvest)
791	fuelwood_biomass(:,:) = zero
792	timber_biomass(:,:) = zero
793	exported_biomass(:,:) = zero
794
795	IF (fuelwood_diameter(ivm) .LT. min_stomate) THEN
796
797	! Grasses and crops
798	IF(harvest_cut(ipts,ivm).EQ.icut_lcc_wood .OR. &
799	harvest_cut(ipts,ivm).EQ.icut_lcc_res) THEN
800
801	! All biomass exported following lcc
802	exported_biomass(:,ilcc) = SUM(harvest_pool_acc(ipts,ivm,:,:),1)
803
804	ELSE
805
806	! All biomass exported following some sort of harvest
807	exported_biomass(:,iharvest) = SUM(harvest_pool_acc(ipts,ivm,:,:),1)
808
809	ENDIF
810
811	! Attribute the harvested biomass to pools with a different
812	! turnover time. The turnover times are PFT depend which
813	! accounts for the different uses different PFT's are put to.
814	! Grasses and crops do not provide fuelwood. They do have
815	! a very short turnover time. Medium and long-lived are
816	! nevertheless considered in case grasses and crops become
817	! to be used in longer-lived products (e.g. bio-plastics)
818	! It also makes the code more consistent. The values stored in
819	! harvest_pool_acc are absolute numbers gC (for that pixel), hence,
820	! the spatial extent of the LCC change, harvest or thinning has
821	! already been accounted for.
822	prod_s(ipts,0,:,:,ilct) = prod_s(ipts,0,:,:,ilct) + &
823	(coeff_lcchange_s(ivm) * exported_biomass(:,:))
824	prod_m(ipts,0,:,:,ilct) = prod_m(ipts,0,:,:,ilct) + &
825	(coeff_lcchange_m(ivm) * exported_biomass(:,:))
826	prod_l(ipts,0,:,:,ilct) = prod_l(ipts,0,:,:,ilct) + &
827	(coeff_lcchange_l(ivm) * exported_biomass(:,:))
828
829	! Calculate how much C and N is released back into the
830	! atmosphere the same year (needed for AR6 output)
831	! Note: when nshort=1, this makes completely sense. For nshort>1,
832	! this only accounts for the carbon release happening within the
833	! the year of harvest.
834	flux_s_pft(ipts,ivm,:,ilct) = SUM(harvest_pool_acc(ipts,ivm,:,:),dim=1) * &
835	coeff_lcchange_s(ivm) / nshort
836
837	ELSE
838
839	! Forests
840	! This is a woody PFT, check the diameter of the harvested wood
841	IF(harvest_cut(ipts,ivm).EQ.icut_lcc_wood .OR. &
842	harvest_cut(ipts,ivm).EQ.icut_lcc_res) THEN
843
844	! All biomass exported following lcc
845	timber_biomass(:,ilcc) = SUM(harvest_pool_acc(ipts,ivm,:,:),1)
846
847	! Use the dimensions of the stems to move the biomass into
848	! the right product pool
849	DO idia = 1,ndia_harvest
850
851	IF (harvest_pool_bound(idia) .LT. fuelwood_diameter(ivm)) THEN
852
853	! If the diameter of the harvest pool is below the treshold
854	! the wood is used as fuelwood and will end up in the short
855	! lived pool
856	fuelwood_biomass(:,ilcc) = fuelwood_biomass(:,ilcc) + &
857	harvest_pool_acc(ipts,ivm,idia,:)
858
859	ENDIF
860
861	ENDDO ! harvest diameter
862
863	ELSE
864
865	! All biomass exported following some sort of harvest
866	timber_biomass(:,iharvest) = SUM(harvest_pool_acc(ipts,ivm,:,:),1)
867
868	! Use the dimensions of the stems to move the biomass into
869	! the right product pool
870	DO idia = 1,ndia_harvest
871
872	IF (harvest_pool_bound(idia) .LT. fuelwood_diameter(ivm)) THEN
873
874	! If the diameter of the harvest pool is below the treshold
875	! the wood is used as fuelwood and will end up in the short
876	! lived pool
877	fuelwood_biomass(:,iharvest) = fuelwood_biomass(:,iharvest) + &
878	harvest_pool_acc(ipts,ivm,idia,:)
879
880	ENDIF
881
882	ENDDO ! harvest diameter
883
884	ENDIF
885
886	! Calculate the long lived pools as the residual of all
887	! wood mines the fuel wood.
888	timber_biomass(:,:) = timber_biomass(:,:) - fuelwood_biomass(:,:)
889
890	! Attribute the harvested biomass to pools with a different
891	! turnover time. The turnover times are PFT depend which
892	! accounts for the different uses different PFT's are put to i.e.
893	! wood from forests vs grass or crops. The values stored in
894	! harvest_pool_acc are absolute numbers gC (for that pixel), hence,
895	! the spatial extent of the LCC change, harvest or thinning has
896	! already been accounted for.
897	prod_s(ipts,0,:,:,ilct) = prod_s(ipts,0,:,:,ilct) + fuelwood_biomass(:,:)
898
899	! Calculate how much C and N is released back into the
900	! atmosphere the same year (needed for AR6 output)
901	flux_s_pft(ipts,ivm,:,ilct) = SUM(prod_s(ipts,0,:,:,ilct),dim=2) + &
902	(coeff_lcchange_s(ivm) * SUM(exported_biomass(:,:),dim=2)) / nshort
903
904	IF (coeff_lcchange_m(ivm) + coeff_lcchange_l(ivm) .GT. &
905	min_stomate) THEN
906
907	! If the condition is satisfied, calculate the fractions
908	! if it is not satisfied there is no medium and long lived
909	! pool for the products. Nothing should then be done. Avoid
910	! divide by zero
911	prod_m(ipts,0,:,:,ilct) = prod_m(ipts,0,:,:,ilct) + &
912	(coeff_lcchange_m(ivm)/ &
913	(coeff_lcchange_m(ivm)+coeff_lcchange_l(ivm)) * &
914	timber_biomass(:,:))
915	prod_l(ipts,0,:,:,ilct) = prod_l(ipts,0,:,:,ilct) + &
916	(coeff_lcchange_l(ivm)/ &
917	(coeff_lcchange_m(ivm)+coeff_lcchange_l(ivm)) * &
918	timber_biomass(:,:))
919
920	ENDIF ! Check coeff_lcchange
921
922	ENDIF ! grasses/crops vs forests
923
924	ENDDO ! loop over PFTs
925
926	DO ilct = 1,nlctypes
927
928	! Check whether all the harvest was used as is to be expected
929	! The precision we are aiming for is 10-8 per m-2. The harvest pool
930	! is expressed for the whole pixel
931	residual(:) = harvest_acc_tmp(:,ilct) - &
932	SUM(prod_s(ipts,0,:,:,ilct),2) - &
933	SUM(prod_m(ipts,0,:,:,ilct),2) - &
934	SUM(prod_l(ipts,0,:,:,ilct),2)
935
936	DO iele = 1,nelements
937
938	IF (ABS(residual(iele)/area(ipts)) .GT. min_stomate) THEN
939
940	WRITE (numout,*) 'Error: some harvest in '//&
941	'sapiens_product_use has not been attributed'
942	WRITE (numout,*) 'Residual in pixel, iele, ', &
943	iele, ipts, residual
944	CALL ipslerr_p (3,'sapiens_product_use', &
945	'basic_product_use','some harvest has not been attributed','')
946
947	ELSE
948
949	! We could still have residual within the precision at the
950	! m2 level but outside the precision at the pixel level.
951	! Move the residual to the short lived pool
952	IF ( prod_s(ipts,0,iele,iharvest,ilct).GT.zero) THEN
953
954	! If there was harvest put the residual in the product pool of
955	! wood harvest
956	prod_s(ipts,0,iele,iharvest,ilct) = prod_s(ipts,0,iele,iharvest,ilct) + &
957	residual(iele)
958
959	ELSE
960
961	! There was no harvest but there is some residual so the residual
962	! should come from land cover change
963	prod_s(ipts,0,iele,ilcc,ilct) = prod_s(ipts,0,iele,ilcc,ilct) + &
964	residual(iele)
965
966	ENDIF
967
968	ENDIF
969
970	END DO ! iele
971
972	END DO ! land cover types
973
974	! All the carbon in the harvest pool was allocated to the product
975	! pools the harvest pool can be set to zero so we can us it with
976	! confidence in the mass balance calculation
977	harvest_pool_acc(ipts,:,:,iele) = zero
978	harvest_type(ipts,:) = zero
979	harvest_cut(ipts,:) = zero
980	harvest_area(ipts,:) = zero
981
982	!! 2.3 Update the long-turnover pool content
983	! Update the long-turnover pool content following flux emission
984	! of a linear decay (1/longivety) of the initial carbon input.
985	! This must be implemented with a counter going from ::nlong to
986	! zero because this way it takes ::nlong years before the intial
987	! pool has been entirely consumed. If a more intuitive approach
988	! would have been used in which the counter goes from zero to
989	! ::nlong, the pool would be depleted in the first year.
990	! Update the long-turnover pool content following flux emission.
991	DO iage = nlong,2,-1
992
993	! below codes are a simple bookkeeping one, note that flux_l contains
994	! the constant annual value that needs to be removed every year.
995	! The key is therefore to just shift their positions
996
997	! Suppose nlong = 10; logic of the following lines are:
998	!
999	! flux_prod = flux_prod + flux(10)
1000	! prod(10) = prod(9) - flux(9) !shift the 9th to 10th
1001	! flux(10) = flux(9) ! this is to prepare for next year
1002
1003	! flux_prod = flux_prod + flux(9)
1004	! prod(9) = prod(8) - flux(8) !shift the 8th to 9th
1005	! flux(9) = flux(8)
1006
1007	! ...
1008
1009	! flux_prod = flux_prod + flux(2)
1010	! prod(2) = prod(1) - flux(1) !shift the 1st to 2nd
1011	! flux(2) = flux(1)
1012
1013	flux_prod_l(ipts,:,:,:) = flux_prod_l(ipts,:,:,:) + flux_l(ipts,iage,:,:,:)
1014	prod_l(ipts,iage,:,:,:) = prod_l(ipts,iage-1,:,:,:) - flux_l(ipts,iage-1,:,:,:)
1015	flux_l(ipts,iage,:,:,:) = flux_l(ipts,iage-1,:,:,:)
1016
1017	ENDDO
1018
1019	! Treat the first year separately
1020	flux_prod_l(ipts,:,:,:) = flux_prod_l(ipts,:,:,:) + flux_l(ipts,1,:,:,:)
1021
1022	! Each year 1/::nlong of the initial pool is decomposed. This fixed
1023	! amount of decomposition is moved to a higher age class each year.
1024	! The units of ::flux_m are gC (absolute number as is the product
1025	! pool). At this point in the calculation, the flux is NOT a flux
1026	! yet but is also a pool. Note harvest pool generated this year only
1027	! starts to decompose the next year.
1028	flux_l(ipts,1,:,:,:) = (un/nlong) * prod_l(ipts,0,:,:,:)
1029	prod_l(ipts,1,:,:,:) = prod_l(ipts,0,:,:,:)
1030	prod_l(ipts,0,:,:,:) = zero
1031
1032	!! 2.4 Update the medium-turnover pool content
1033	! Update the medium-turnover pool content following flux emission.
1034	! For comments see the section above. The concept is identical.
1035	DO iage = nmedium,2,-1
1036	flux_prod_m(ipts,:,:,:) = flux_prod_m(ipts,:,:,:) + flux_m(ipts,iage,:,:,:)
1037	prod_m(ipts,iage,:,:,:) = prod_m(ipts,iage-1,:,:,:) - flux_m(ipts,iage-1,:,:,:)
1038	flux_m(ipts,iage,:,:,:) = flux_m(ipts,iage-1,:,:,:)
1039	ENDDO
1040
1041	flux_prod_m(ipts,:,:,:) = flux_prod_m(ipts,:,:,:) + flux_m(ipts,1,:,:,:)
1042	flux_m(ipts,1,:,:,:) = (un/nmedium) * prod_m(ipts,0,:,:,:)
1043	prod_m(ipts,1,:,:,:) = prod_m(ipts,0,:,:,:)
1044	prod_m(ipts,0,:,:,:) = zero
1045
1046	!! 2.5 Update the short-turnover pool
1047	! The length of the short-turnover pool will most likely be
1048	! less than length of the regular loop. Different from long- and
1049	! medium-residence pool, we suppose decomposition of short pool
1050	! start from the year of harvest, rather than the next year of
1051	! harvest year.
1052	IF (nshort .GE. 3) THEN
1053
1054	! note we start from nshort-1 here,because starting from the harvest
1055	! year, it takes 'nshort-1' years for the pool to completely return to
1056	! the atmosphere.
1057	DO iage = nshort-1,2,-1
1058	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + flux_s(ipts,iage,:,:,:)
1059	prod_s(ipts,iage,:,:,:) = prod_s(ipts,iage-1,:,:,:) - flux_s(ipts,iage-1,:,:,:)
1060	flux_s(ipts,iage,:,:,:) = flux_s(ipts,iage-1,:,:,:)
1061	ENDDO
1062
1063	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + flux_s(ipts,1,:,:,:)
1064	flux_s(ipts,1,:,:,:) = (un/nshort) * prod_s(ipts,0,:,:,:)
1065	! this line accounts for the decomposition happening at the year
1066	! of harvest.
1067	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + flux_s(ipts,1,:,:,:)
1068	prod_s(ipts,1,:,:,:) = prod_s(ipts,0,:,:,:) - flux_s(ipts,1,:,:,:)
1069	prod_s(ipts,0,:,:,:) = zero
1070
1071	ELSEIF (nshort .EQ. 2) THEN
1072	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + flux_s(ipts,1,:,:,:)
1073	flux_s(ipts,1,:,:,:) = (un/nshort) * prod_s(ipts,0,:,:,:)
1074	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + flux_s(ipts,1,:,:,:)
1075	prod_s(ipts,1,:,:,:) = prod_s(ipts,0,:,:,:) - flux_s(ipts,1,:,:,:)
1076	prod_s(ipts,0,:,:,:) = zero
1077
1078	ELSEIF (nshort .EQ. 1) THEN
1079
1080	! Everything decomposes in a single year. The short product
1081	! pool remains empty.
1082	flux_prod_s(ipts,:,:,:) = flux_prod_s(ipts,:,:,:) + prod_s(ipts,0,:,:,:)
1083	prod_s(ipts,0,:,:,:) = zero
1084
1085	ELSE
1086
1087	WRITE(numout,*) 'Error: flaw in the logic of IF-construct '//&
1088	'(basic_product_use in sapiens_product_use.f90)'
1089	CALL ipslerr_p(3,'ERROR: product_use',&
1090	'Logic flaw in if-construct','','')
1091
1092	ENDIF ! number of ages in short-turnover product pool
1093
1094	ENDDO ! #pixels - spatial domain
1095
1096
1097	!! 3. Check numerical consistency of this routine
1098
1099	IF (err_act.GT.1) THEN
1100
1101	! 3.1 Mass balance closure
1102	! 3.1.1 Calculate final biomass
1103	! The biomass harvest pool shouldn't be multiplied by veget_max
1104	! its in gC pixel-1 so divide by area to get a value in
1105	! gC m-2. Account for the harvest and wood product pools. The harvest
1106	! pool should be empty, it was added for completeness.
1107	pool_end = zero
1108	DO iele=1,nelements
1109	pool_end(:,1,iele) = &
1110	( SUM(SUM(harvest_pool_acc(:,:,:,iele),3),2) + &
1111	SUM(SUM(SUM(prod_l(:,:,iele,:,:),2),2),2) + &
1112	SUM(SUM(SUM(prod_m(:,:,iele,:,:),2),2),2) + &
1113	SUM(SUM(SUM(prod_s(:,:,iele,:,:),2),2),2) ) / area(:)
1114	ENDDO
1115
1116	!! 3.1.2 Calculate mass balance
1117	! Common processes
1118	DO iele=1,nelements
1119	check_intern(:,1,iland2atm,iele) = &
1120	-un * SUM(SUM(flux_prod_s(:,iele,:,:) + &
1121	flux_prod_m(:,iele,:,:) + &
1122	flux_prod_l(:,iele,:,:),2),2 ) / &
1123	area(:)
1124	check_intern(:,1,ipoolchange,iele) = &
1125	-un * (pool_end(:,1,iele) - &
1126	pool_start(:,1,iele))
1127	ENDDO
1128
1129	closure_intern = zero
1130	DO imbc = 1,nmbcomp
1131	DO iele=1,nelements
1132	closure_intern(:,1,iele) = closure_intern(:,1,iele) + &
1133	check_intern(:,1,imbc,iele)
1134	ENDDO
1135	ENDDO
1136
1137	! 3.1.3 Check mass balance closure
1138	CALL check_mass_balance("basic_product_use", closure_intern, npts, &
1139	pool_end, pool_start, veget_max, 'products')
1140
1141	ENDIF ! err_act.GT.1
1142
1143	!! 4. Convert the pools into fluxes and aggregate
1144
1145	! Convert pools into fluxes
1146	! The pools are given in the absolute amount of carbon for the
1147	! whole pixel and the harvest was accumulated over the year.
1148	! The unit is thus gc pixel-1 year-1. These fluxes should be
1149	! expressed as gC m-2 dt-1. Therefore, divide by the area
1150	! of the pixel and the number of time steps within a year
1151	DO iele = 1,nelements
1152
1153	DO ilan = 1,nlanduse
1154
1155	DO ilct = 1,nlctypes
1156
1157	flux_prod_s(:,iele,ilan,ilct) = flux_prod_s(:,iele,ilan,ilct) / area(:)
1158	flux_prod_m(:,iele,ilan,ilct) = flux_prod_m(:,iele,ilan,ilct) / area(:)
1159	flux_prod_l(:,iele,ilan,ilct) = flux_prod_l(:,iele,ilan,ilct) / area(:)
1160
1161	! Calculate aggregate values
1162	! flux_prod_total is now in gC m-2 dt-1. The prod_x_total
1163	! are in gC m-2
1164	flux_prod_total(:,iele,ilan,ilct) = flux_prod_s(:,iele,ilan,ilct) + &
1165	flux_prod_m(:,iele,ilan,ilct) + flux_prod_l(:,iele,ilan,ilct)
1166	prod_s_total(:,iele,ilan,ilct) = SUM(prod_s(:,:,iele,ilan,ilct),2)
1167	prod_m_total(:,iele,ilan,ilct) = SUM(prod_m(:,:,iele,ilan,ilct),2)
1168	prod_l_total(:,iele,ilan,ilct) = SUM(prod_l(:,:,iele,ilan,ilct),2)
1169
1170	END DO
1171
1172	END DO
1173
1174	! AR6 output
1175	flux_s_pft(:,:,iele,:) = flux_s_pft(:,:,iele,:) / one_year * dt_days
1176
1177	END DO
1178
1179	IF (printlev.GE.4) WRITE(numout,*) 'Leaving dimensional product use'
1180
1181
1182	END SUBROUTINE dim_product_use
1183
1184	!! ================================================================================================================================
1185	!! SUBROUTINE : product_init
1186	!!
1187	!>\BRIEF Initialize the product use variables when product use is not
1188	!! calculated
1189	!!
1190	!! DESCRIPTION : Initialize the product use variables when product use is not
1191	!! calculated
1192	!!
1193	!! RECENT CHANGE(S) : None
1194	!!
1195	!! MAIN OUTPUT VARIABLE(S) : prod_s, prod_m, prod_l, flux_s, flux_m, flux_l,
1196	!! flux_prod_s flux_prod_m and flux_prod_l
1197	!!
1198	!! REFERENCES : None
1199	!!
1200	!! FLOWCHART : None
1201	!!
1202	!_ ================================================================================================================================
1203
1204	SUBROUTINE product_init(flux_prod_s, flux_prod_m, flux_prod_l, &
1205	prod_s_total, prod_m_total, prod_l_total, flux_s, &
1206	flux_m, flux_l, flux_prod_total, flux_s_pft, &
1207	prod_s, prod_m, prod_l)
1208
1209	IMPLICIT NONE
1210
1211
1212	!! 0. Variable and parameter declaration
1213
1214	!! 0.1 Input variables
1215
1216	!! 0.2 Output variables
1217
1218	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_s !! C-released during first years (short term)
1219	!! following land cover change
1220	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
1221	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_m !! Total annual release from decomposition of
1222	!! the medium-lived product pool
1223	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
1224	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_l !! Total annual release from decomposition of
1225	!! the long-lived product pool
1226	!! @tex ($gC m^{-1} dt^{-1}$) @endtex
1227	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: prod_s_total !! Total products remaining in the short-lived
1228	!! pool after the annual decomposition
1229	!! @tex $(gC)$ @endtex
1230	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: prod_m_total !! Total products remaining in the medium-lived
1231	!! pool after the annual decomposition
1232	!! @tex $(gC)$ @endtex
1233	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: prod_l_total !! Total products remaining in the long-lived
1234	!! pool after the annual release
1235	!! @tex $(gC)$ @endtex
1236	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_prod_total !! Total flux from decomposition of the short,
1237	!! medium and long lived product pools
1238	!! @tex $(gC m^{-1} dt^{-1})$ @endtex
1239	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: flux_s_pft !! Total flux from decomposition the same year of the lcc and harvest
1240	!! @tex $(gC m^{-1} dt^{-1})$ @endtex
1241
1242	!! 0.3 Modified variables
1243
1244	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: flux_s !! Annual release from the short-lived product
1245	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
1246	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: flux_m !! Annual release from the medium-lived product
1247	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
1248	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: flux_l !! Annual release from the long-lived product
1249	!! pool @tex $(gC) pixel^{-1} year^{-1}$ @endtex
1250	REAL(r_std), DIMENSION(:,0:,:,:,:), INTENT(inout) :: prod_s !! Short-lived product pool after the annual
1251	!! release of each compartment (short + 1 :
1252	!! input from year of land cover change)
1253	!! @tex ($gC$) @endtex
1254	REAL(r_std), DIMENSION(:,0:,:,:,:), INTENT(inout) :: prod_m !! Medium-lived product pool after the annual
1255	!! release of each compartment (medium + 1 :
1256	!! input from year of land cover change)
1257	!! @tex ($gC$) @endtex
1258	REAL(r_std), DIMENSION(:,0:,:,:,:), INTENT(inout) :: prod_l !! long-lived product pool after the annual
1259	!! release of each compartment (long + 1 :
1260	!! input from year of land cover change)
1261	!! @tex ($gC$) @endtex
1262
1263	!! 0.4 Local variables
1264
1265	!_ ================================================================================================================================
1266
1267	IF (printlev.GE.2) WRITE(numout,*) 'Entering product init'
1268
1269	!! 1. initialization
1270	prod_s(:,0,:,:,:) = zero
1271	prod_m(:,0,:,:,:) = zero
1272	prod_l(:,0,:,:,:) = zero
1273	flux_prod_s(:,:,:,:) = zero
1274	flux_prod_m(:,:,:,:) = zero
1275	flux_prod_l(:,:,:,:) = zero
1276	prod_s_total(:,:,:,:) = zero
1277	prod_m_total(:,:,:,:) = zero
1278	prod_l_total(:,:,:,:) = zero
1279	flux_s_pft(:,:,:,:) = zero
1280	flux_prod_total(:,:,:,:) = zero
1281
1282	IF (printlev.GE.4) WRITE(numout,*) 'Leaving product init'
1283
1284	END SUBROUTINE product_init
1285
1286	END MODULE sapiens_product_use
1287
1288
1289

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

Download in other formats: