Context Navigation

stomate_npp.f90 @ 7346

Last change on this file since 7346 was 5488, checked in by chunjing.qiu, 6 years ago
C balance checked
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 39.0 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_npp
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 This modules calculates NPP: Maintenance and growth respiration
10	!!
11	!!\n DESCRIPTION: We calculate first the maintenance respiration. This is substracted from the
12	!! allocatable biomass (and from the present biomass if the GPP is too low).\n
13	!! Of the rest, a part is lost as growth respiration, while the other part is
14	!! effectively allocated.
15	!!
16	!! RECENT CHANGE(S): None
17	!!
18	!! REFERENCE(S) :
19	!!
20	!! SVN :
21	!! $HeadURL$
22	!! $Date$
23	!! $Revision$
24	!! \n
25	!_ ================================================================================================================================
26
27	MODULE stomate_npp
28
29	! modules used:
30	USE xios_orchidee
31	USE ioipsl_para
32	USE stomate_data
33	USE constantes
34	USE constantes_soil
35	USE pft_parameters
36	USE crop_alloc
37
38	IMPLICIT NONE
39
40	! private & public routines
41
42	PRIVATE
43	PUBLIC npp_calc,npp_calc_clear
44
45	LOGICAL, SAVE :: firstcall_npp = .TRUE. !! first call
46	!$OMP THREADPRIVATE(firstcall_npp)
47
48	CONTAINS
49
50	!! ================================================================================================================================
51	!! SUBROUTINE : npp_calc_clear
52	!!
53	!>\BRIEF : Set the flag ::firstcall_npp to .TRUE. and as such activate section
54	!! 1.1 of the subroutine npp_calc (see below).\n
55	!_ ================================================================================================================================
56
57	SUBROUTINE npp_calc_clear
58	firstcall_npp=.TRUE.
59	END SUBROUTINE npp_calc_clear
60
61
62
63
64
65	!! ================================================================================================================================
66	!! SUBROUTINE : npp_calc
67	!!
68	!>\BRIEF Calculate NPP as the difference between GPP and respiration (= growth + maintenance respiration).
69	!! Update biomass of all compartments after calculating respiration and allocation.
70	!!
71	!!
72	!! DESCRIPTION : NPP is calculated from three components: Gross Primary Productivity (GPP), maintenance respiration
73	!! and growth respiration (all in @tex $ gC.m^{-2}dt^{-1} $ @endtex), following the convention that positive fluxes denote
74	!! fluxes plants to the atmosphere. GPP is the input variable from which, in the end, NPP or total allocatable biomass
75	!! @tex $(gC.m^{-2}dt^{-1}))$ @endtex is calculated. Net primary production is then calculated as:\n
76	!! NPP = GPP - growth_resp - maint-resp [eq. 1]\n
77	!!
78	!! The calculation of maintenance respiration is done in routine stomate_resp.f90. Maintenance respiration is calculated for
79	!! the whole plant and is therefore removed from the total allocatable biomass. In order to prevent all allocatable biomass
80	!! from being used for maintenance respiration, a limit fraction of total allocatable biomass, tax_max, is defined (in
81	!! variables declaration). If maintenance respiration exceeds tax_max (::bm_tax_max), the maximum allowed allocatable biomass
82	!! will be respired and the remaining respiration, required in excess of tax_max, is taken out from tissues already present in
83	!! the plant (biomass).\n
84	!!
85	!! After total allocatable biomass has been updated by removing maintenance respiration, total allocatable biomass is distributed
86	!! to all plant compartments according to the f_alloc fractions calculated in stomate_alloc.f90.\n
87	!!
88	!! Growth respiration is calculated as a fraction of allocatable biomass for each part of the plant. The fraction coefficient
89	!! ::frac_growth_resp is defined in stomate_constants.f90 and is currently set to be the same for all plant compartments.
90	!! Allocatable biomass of all plant compartments are updated by removing what is lost through growth respiration. Net allocatable
91	!! biomass (total allocatable biomass after maintenance and growth respiration) is added to the current biomass for each plant
92	!! compartment.
93	!!
94	!! Finally, leaf age and plant age are updated. Leaf age is described with the concept of "leaf age classes". A number of leaf
95	!! age classes (nleafages) is defined in stomate_constants.f90. Each leaf age class contains a fraction (::leaf_frac) of the
96	!! total leaf biomass. When new biomass is added to leaves, the age of the biomass in the youngest leaf age class is decreased.
97	!! The fractions of leaves in the other leaf ages classes are also updated as the total biomass has increased. Plant age is
98	!! updated first by increasing the age of the previous biomass by one time step, and then by adjusting this age as the average
99	!! of the ages of the previous and the new biomass.
100	!!
101	!! RECENT CHANGE(S): None
102	!!
103	!! MAIN OUTPUT VARIABLE(S): ::npp
104	!!
105	!! REFERENCE(S) :
106	!! - F.W.T.Penning De Vries, A.H.M. Brunsting, H.H. Van Laar. 1974. Products, requirements and efficiency of biosynthesis a
107	!! quantitative approach. Journal of Theoretical Biology, Volume 45, Issue 2, June 1974, Pages 339-377.
108	!!
109	!! FLOWCHART :
110	!! \latexonly
111	!! \includegraphics[scale=0.14]{stomate_npp_flow.jpg}
112	!! \endlatexonly
113	!! \n
114	!_ ================================================================================================================================
115
116	SUBROUTINE npp_calc (npts, dt, &
117	PFTpresent, &
118	t2m, tsoil, lai, rprof, &
119	gpp, f_alloc, bm_alloc, resp_maint_part,&
120	biomass, leaf_age, leaf_frac, age, &
121	resp_maint, resp_growth, npp, &
122	!!!! crop variables
123	! for crop bm_alloc
124	in_cycle, deltai, dltaisen, ssla, pgrain, deltgrain, reprac, &
125	nger, nlev, ndrp, nlax, nmat, nrec, &
126	c_reserve, c_leafb, slai, tday_counter, veget_max, &
127	!!!! end crop variables, xuhui
128	!gmjc
129	sla_calc, sla_age1,N_limfert)
130	!end gmjc
131	!! 0 Variable and parameter declaration
132
133	!! 0.1 Input variables
134
135	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of pixels (unitless)
136	REAL(r_std), INTENT(in) :: dt !! Time step (days)
137	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: PFTpresent !! PFT exists (true/false)
138	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m !! Temperature at 2 meter (K)
139	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: tsoil !! Soil temperature of each soil layer (K)
140	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lai !! PFT leaf area index (unitless)
141	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: rprof !! PFT root depth as calculated in stomate.f90
142	!! from root profile parameter humcste (m)
143	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: gpp !! PFT gross primary productivity
144	!! @tex $(gC.m^{-2}dt^{-1})$ @endtex
145	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: f_alloc !! Fraction of total allocatable biomass that
146	!! goes into each plant part (unitless)
147	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: resp_maint_part !! Maintenance respiration of different plant
148	!! parts @tex $(gC.m^{-2}dt^{-1})$ @endtex
149	!!!!! crop var
150	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: in_cycle
151	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: deltai
152	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: dltaisen
153	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ssla
154	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: pgrain
155	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: deltgrain
156	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: reprac
157	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nger
158	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nlev
159	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: ndrp
160	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nlax
161	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nmat
162	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nrec
163	INTEGER(i_std), INTENT(in) :: tday_counter
164	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max
165	!!!!! xuhui
166
167	!! 0.2 Output variables
168
169	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_maint !! PFT maintenance respiration
170	!! @tex $(gC.m^{-2}dt^{-1})$ @endtex
171	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_growth !! PFT growth respiration
172	!! @tex $(gC.m^{-2}dt^{-1})$ @endtex
173	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: npp !! PFT net primary productivity
174	!! @tex $(gC.m^{-2}dt^{-1})$ @endtex
175	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out) :: bm_alloc !! PFT biomass increase, i.e. NPP per plant part
176	!! @tex $(gC.m^{-2}dt^{-1})$ @endtex
177
178	!! 0.3 Modified variables
179
180	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass !! PFT total biomass of each plant part
181	!! @tex $(gC.m^{-2})$ @endtex
182	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_age !! PFT age of different leaf age classes (days)
183	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac !! PFT fraction of total leaves in leaf age
184	!! class (unitless)
185	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age !! PFT age (years)
186	!gmjc
187	! N fertilization limitation factor for grassland Vcmax and SLA
188	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: N_limfert
189	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: sla_calc
190	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: sla_age1
191	!end gmjc
192	!!!!! crop var
193	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: c_reserve
194	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: c_leafb
195	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: slai
196	!!!!! xuhui
197	!! 0.4 Local variables
198	!gmjc
199	! weighted leaf age (leaf_age * fraction of the age class)
200	REAL(r_std), DIMENSION(npts,nvm) :: leaf_age_w
201	REAL(r_std), DIMENSION(npts,nvm) :: sla_age2
202	REAL(r_std), DIMENSION(npts,nvm) :: sla_age3
203	REAL(r_std), DIMENSION(npts,nvm) :: sla_age4
204	! SLA max and SLA min affected by N fertilization
205	REAL(r_std), DIMENSION(npts,nvm) :: sla_max_Nfert
206	REAL(r_std), DIMENSION(npts,nvm) :: sla_min_Nfert
207	!end gmjc
208	!BEGINNVADD
209	! npp above (gC/m**2) for the whole plant
210	REAL(r_std), DIMENSION(npts,nvm) :: npp_above
211	! npp below (gC/m**2) for the whole plant
212	REAL(r_std), DIMENSION(npts,nvm) :: npp_below
213	! NPP per plant part
214	REAL(r_std), DIMENSION(npts,nvm,nparts) :: npp_part
215	!ENDNVADD
216	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: z_soil !! Soil levels representing soil depth (m)
217	!$OMP THREADPRIVATE(z_soil)
218	REAL(r_std), DIMENSION(npts,nvm) :: t_root !! Root temperature (convolution of root and
219	!! soil temperature profiles)(K)
220	REAL(r_std), DIMENSION(npts,nvm,nparts) :: coeff_maint !! PFT maintenance respiration coefficients of
221	!! different plant compartments at 0 deg C
222	!! @tex $(g.g^{-1}dt^{-1})$ @endtex
223	REAL(r_std), DIMENSION(npts,nparts) :: t_maint !! Temperature which is pertinent for maintenance
224	!! respiration, which is air/root temperature for
225	!! above/below-ground compartments (K)
226	REAL(r_std), DIMENSION(npts) :: rpc !! Scaling factor for integrating vertical soil
227	!! profiles (unitless)
228	REAL(r_std), DIMENSION(npts) :: tl !! Long term annual mean temperature (C)
229	REAL(r_std), DIMENSION(npts) :: slope !! Slope of maintenance respiration coefficient
230	!! (1/K)
231	REAL(r_std), DIMENSION(npts,nparts) :: resp_growth_part !! Growth respiration of different plant parts
232	!! @tex $(gC.m^{-2}dt^{-1})$ @endtex
233	REAL(r_std), DIMENSION(npts,nvm) :: bm_alloc_tot !! Allocatable biomass for the whole plant
234	!! @tex $(gC.m^{-2})$ @endtex
235	REAL(r_std), DIMENSION(npts) :: bm_add !! Biomass increase @tex $(gC.m^{-2})$ @endtex
236	REAL(r_std), DIMENSION(npts) :: bm_new !! New biomass @tex $(gC.m^{-2})$ @endtex
237	REAL(r_std), DIMENSION(npts,nvm) :: leaf_mass_young !! Leaf mass in youngest age class
238	!! @tex $(gC.m^{-2})$ @endtex
239	REAL(r_std), DIMENSION(npts,nvm) :: lm_old !! Leaf mass after maintenance respiration
240	!! @tex $(gC.m^{-2})$ @endtex
241	REAL(r_std), DIMENSION(npts,nvm) :: bm_create !! Biomass created when biomass<0 because of dark
242	!! respiration @tex $(gC.m^{-2})$ @endtex
243	REAL(r_std), DIMENSION(npts) :: bm_tax_max !! Maximum part of allocatable biomass used for
244	!! respiration @tex $(gC.m^{-2})$ @endtex
245	REAL(r_std), DIMENSION(npts) :: bm_pump !! Biomass that remains to be taken away
246	!! @tex $(gC.m^{-2})$ @endtex
247	INTEGER(i_std) :: i,j,k,l,m !! Indeces(unitless)
248	INTEGER(i_std) :: ier !! Error handling
249
250	!_ ================================================================================================================================
251
252	IF (printlev>=3) WRITE(numout,*) 'Entering npp'
253
254	!! 1. Initializations
255
256	!! 1.1 First call
257	IF ( firstcall_npp ) THEN
258
259	!! 1.1.1 Soil levels
260	! Get the depth of the different soil layers (number of layers=nbdl)
261	! previously calculated as variable diaglev in routines sechiba.f90 and slowproc.f90
262	ALLOCATE(z_soil(0:nslm), stat=ier)
263	IF ( ier /= 0 ) CALL ipslerr_p(3,'npp_calc','Pb in allocate of z_soil','','')
264
265	z_soil(0) = zero
266	z_soil(1:nslm) = diaglev(1:nslm)
267
268	!! 1.1.2 Output message
269	! Write message including value used for tax_max
270	WRITE(numout,*) 'npp:'
271
272	WRITE(numout,*) ' > max. fraction of allocatable biomass used for'// &
273	' maint. resp.:', tax_max
274
275	firstcall_npp = .FALSE.
276
277	ENDIF ! End if first call
278
279	!! 1.2 Set output variables to zero
280	bm_alloc(:,:,:,:) = zero
281	resp_maint(:,:) = zero
282	resp_growth(:,:) = zero
283	npp(:,:) = zero
284	!BEGINNVADD
285	npp_above(:,:) = zero
286	npp_below(:,:) = zero
287	npp_part(:,:,:) = zero
288	!ENDNVADD
289
290	!! 1.3 Total allocatable biomass
291	! total allocatable biomass during this time step determined from GPP.
292	! GPP was calculated as CO2 assimilation in enerbil.f90
293	!WRITE(numout,*) 'zd bm_alloc 1','bm_alloc_tot(1,10)',bm_alloc_tot(1,10),'gpp(1,10)',gpp(1,10),'dt',dt
294	bm_alloc_tot(:,:) = gpp(:,:) * dt
295	!WRITE(numout,*) 'zd bm_alloc 2','bm_alloc_tot(1,10)',bm_alloc_tot(1,10)
296
297
298	! WRITE(numout,*) 'biomass(1,12:14,:,icarbon) before PFT loop:',biomass(1,12:14,:,icarbon)
299	!! 3. Calculate maintenance and growth respiration
300	! First, total maintenance respiration for the whole plant is calculated by summing maintenance
301	! respiration of the different plant compartments. Then, maintenance respiration is subtracted
302	! from whole-plant allocatable biomass (up to a maximum fraction of the total allocatable biomass).
303	! Growth respiration is then calculated for each plant compartment as a fraction of remaining
304	! allocatable biomass for this compartment. NPP is calculated by substracting total autotrophic
305	! respiration from GPP i.e. NPP = GPP - maintenance resp - growth resp.
306	!WRITE(numout,*) 'zd nppcalc1 biomass(1,10,ileaf,icarbon)',biomass(1,10,ileaf,icarbon),'bm_pump(1)',bm_pump(1),'resp_maint_part(1,10,ileaf)',resp_maint_part(1,10,ileaf),'resp_maint(1,10)',resp_maint(1,10)
307	DO j = 2,nvm ! Loop over # of PFTs
308
309	!! 3.1 Maintenance respiration of the different plant parts
310	! Maintenance respiration of the different plant parts is calculated in
311	! stomate_resp.f90 as a function of the plant's temperature,
312	! the long term temperature and plant coefficients
313	! VPP killer:
314	resp_maint(:,j) = zero
315
316	! Following the calculation of hourly maintenance respiration, verify that
317	! the PFT has not been killed after calcul of resp_maint_part in stomate.
318	DO k= 1, nparts
319	WHERE (PFTpresent(:,j))
320	resp_maint(:,j) = resp_maint(:,j) + resp_maint_part(:,j,k)
321	ENDWHERE
322	ENDDO
323
324	!! 3.2 Substract maintenance respiration from allocatable biomass
325	! The total maintenance respiration calculated in 3.2 is substracted from the newly
326	! produced allocatable biomass (bm_alloc_tot). However, ensure that not all allocatable
327	! biomass is removed by setting a maximum to the fraction of allocatable biomass used
328	! for maintenance respiration: tax_max. If the maintenance respiration is larger than
329	! tax_max,the amount tax_max is taken from allocatable biomass, and the remaining of
330	! maintenance respiration is taken from the tissues themselves (biomass). We suppose
331	! that respiration is not dependent on leaf age -> therefore the leaf age structure is
332	! not changed.
333	! The maximum fraction of allocatable biomass used for respiration is defined as tax_max.
334	! The value of tax_max is set in the declarations section (0.4 Local variables) of this
335	! routine
336	bm_tax_max(:) = tax_max * bm_alloc_tot(:,j)
337
338	! IF (j.EQ.10) WRITE(numout,*) 'zd bm_alloc 3','bm_alloc_tot(1,10)',bm_alloc_tot(1,10)
339	DO i = 1, npts ! Loop over # of pixels
340
341	! If there is enough allocatable biomass to cover maintenance respiration,
342	! then biomass associated with maintenance respiration is removed from allocatable biomass
343	IF ( bm_alloc_tot(i,j) .GT. zero ) THEN
344	IF ( ( resp_maint(i,j) * dt ) .LT. bm_tax_max(i) ) THEN
345
346	bm_alloc_tot(i,j) = bm_alloc_tot(i,j) - resp_maint(i,j) * dt
347
348	! If there is not enough allocatable biomass to cover maintenance respiration, the
349	! - maximum allowed allocatable biomass (bm_tax_max) is removed from allocatable biomass.
350	ELSE
351
352	bm_alloc_tot(i,j) = bm_alloc_tot(i,j) - bm_tax_max(i)
353
354	! ::bm_pump is the amount of maintenance respiration that exceeds the maximum allocatable biomass
355	! This amount of biomass still needs to be respired and will be removed from tissues biomass of each
356	! plant compartment
357	bm_pump(i) = resp_maint(i,j) * dt - bm_tax_max(i)
358
359	! The biomass is removed from each plant compartment tissues as the ratio of the maintenance
360	! respiration of the plant compartment to the total maintenance respiration (resp_maint_part/resp_maint)
361	biomass(i,j,ileaf,icarbon) = biomass(i,j,ileaf,icarbon) - &
362	bm_pump(i) * resp_maint_part(i,j,ileaf) / resp_maint(i,j)
363	biomass(i,j,isapabove,icarbon) = biomass(i,j,isapabove,icarbon) - &
364	bm_pump(i) * resp_maint_part(i,j,isapabove) / resp_maint(i,j)
365	biomass(i,j,isapbelow,icarbon) = biomass(i,j,isapbelow,icarbon) - &
366	bm_pump(i) * resp_maint_part(i,j,isapbelow) / resp_maint(i,j)
367	biomass(i,j,iroot,icarbon) = biomass(i,j,iroot,icarbon) - &
368	bm_pump(i) * resp_maint_part(i,j,iroot) / resp_maint(i,j)
369	biomass(i,j,ifruit,icarbon) = biomass(i,j,ifruit,icarbon) - &
370	bm_pump(i) * resp_maint_part(i,j,ifruit) / resp_maint(i,j)
371	biomass(i,j,icarbres,icarbon) = biomass(i,j,icarbres,icarbon) - &
372	bm_pump(i) * resp_maint_part(i,j,icarbres) / resp_maint(i,j)
373	!BEGINNVADD
374	npp_part (i,j,:)= npp_part (i,j,:) - &
375	bm_pump(i) * resp_maint_part(i,j,:) / resp_maint(i,j)
376	!ENDNVADD
377	ENDIF
378	ELSE
379	biomass(i,j,icarbres,icarbon) = biomass(i,j,icarbres,icarbon) + &
380	bm_alloc_tot(i,j) - resp_maint(i,j) * dt
381	bm_alloc_tot(i,j) = 0.
382	ENDIF ! End if there is enough allocatable biomass to cover maintenance respiration
383
384	ENDDO ! Fortran95: WHERE - ELSEWHERE construct
385	! IF (j.EQ.10) WRITE(numout,*) 'zd bm_alloc 4','bm_alloc_tot(1,10)',bm_alloc_tot(1,10)
386
387
388	!! 3.3 Allocate allocatable biomass to different plant compartments.
389	! The amount of allocatable biomass of each compartment is a fraction according f_alloc of total
390	! allocatable biomass (the f_alloc of the different plant parts are calculated in stomate_alloc.f90)
391	! IF (j.EQ.10) WRITE(numout,*) 'zd bm_alloc 5','bm_alloc(1,10,:,icarbon)',bm_alloc(1,10,:,icarbon),'f_alloc(1,10,:)',f_alloc(1,10,:)
392	IF (ok_LAIdev(j)) THEN
393	! WRITE(numout,*) 'slai before npp_alloc: ',j, 'pft, ', slai(1,j)
394	! WRITE(numout,*) 'ssla before npp_alloc: ',j, 'pft, ', ssla(1,j)
395	! WRITE(numout,*) 'biomass(1,j,:,icarbon): ',j, 'pft, ',biomass(1,j,:,icarbon)
396	! WRITE(numout,*) 'bm_alloc_tot(1,j) before npp_alloc: ',j, 'pft, ', bm_alloc_tot(1,j)
397	! WRITE(numout,*) 'bm_alloc(1,j,:,icarbon) before npp_alloc: ',j, 'pft, ', bm_alloc(1,j,:,icarbon)
398	DO i = 1, npts
399	! we call the crop bm allocation subroutine
400	if ( in_cycle(i,j) ) then
401	write(numout,*) '(i,j)', i, j
402	write(numout,*) 'in_cycle(i,j)',in_cycle(i,j)
403	if (veget_max(i,j) .gt. 0.) then
404	! there will be error in crop_alloc if vegetmax for a certain crop is 0
405	call crop_bmalloc(in_cycle(i, j), &
406	deltai(i, j), &
407	dltaisen(i, j), &
408	ssla(i, j), &
409	pgrain(i, j), &
410	deltgrain(i, j), &
411	reprac(i, j), &
412	nger(i, j), &
413	nlev(i,j), &
414	ndrp(i,j), &
415	nlax(i,j), & ! input
416	nmat(i,j), &
417	nrec(i, j), &
418	! f_crop_recycle(i,j), &
419	bm_alloc_tot(i,j), & ! input
420	biomass(i, j, :, icarbon),& ! input
421	c_reserve(i,j), & ! inout
422	c_leafb(i,j), & ! inout
423	bm_alloc(i, j, :, icarbon), & ! inout
424	SP_densitesem(j), &
425	SP_pgrainmaxi(j), &
426	SP_tigefeuil(j), &
427	SP_slamax(j), &
428	slai(i,j), &
429	tday_counter) ! parameter
430	endif ! j==11
431	endif ! in_cycle
432
433	ENDDO ! npts
434	! WRITE(numout,*) 'bm_alloc(1,j,:,icarbon) after npp_alloc: ',j ,'pft, ',bm_alloc(1,j,:,icarbon)
435	! WRITE(numout,*) 'slai after npp_alloc: ',j , 'pft, ',slai(1,j)
436	ELSE ! natural vegetation (is not crop)
437	IF ( (printlev>=5) .AND. (j .EQ. 10) ) THEN
438	WRITE(numout,) 'gpp(:,10)dt', gpp(:,10)*dt
439	WRITE(numout,*) 'bm_alloc_tot(:,10)', bm_alloc_tot(:,10)
440	WRITE(numout,*) 'bm_alloc(:,10,:,icarbon)', bm_alloc(:,10,:,icarbon)
441	ENDIF
442	DO k = 1, nparts
443	bm_alloc(:,j,k,icarbon) = f_alloc(:,j,k) * bm_alloc_tot(:,j)
444	ENDDO
445	ENDIF
446	! IF (j.EQ.10) WRITE(numout,*) 'zd bm_alloc 6','bm_alloc(1,10,:,icarbon)',bm_alloc(1,10,:,icarbon)
447
448
449	!! 3.4 Calculate growth respiration of each plant compartment.
450	! Growth respiration of a plant compartment is a fraction of the allocatable biomass remaining after
451	! maintenance respiration losses have been taken into account. The fraction of allocatable biomass
452	! removed for growth respiration is the same for all plant compartments and is defined by the parameter
453	! frac_growth_resp in stomate_constants.f90. Allocatable biomass ::bm_alloc is updated as a result of
454	! the removal of growth resp.
455
456	!!! xuhui: note that we should exclude those negative bm_alloc induced by crop remobilizations
457	! WRITE(numout,*) 'bm_alloc(1,j,:,icarbon) after growth_resp: ',j ,'pft, ',bm_alloc(1,j,:,icarbon)
458	DO i=1, npts
459	resp_growth_part(i,:) = zero
460	WHERE (bm_alloc(i,j,:,icarbon) .GT. zero)
461	resp_growth_part(i,:) = frac_growthresp(j) * bm_alloc(i,j,:,icarbon) / dt
462	bm_alloc(i,j,:,icarbon) = ( un - frac_growthresp(j) ) * bm_alloc(i,j,:,icarbon)
463	ENDWHERE
464	ENDDO
465	! WRITE(numout,*) 'bm_alloc(1,j,:,icarbon) after growth_resp: ',j ,'pft, ',bm_alloc(1,j,:,icarbon)
466	! IF (j.EQ.10) WRITE(numout,*) 'zd bm_alloc 7','bm_alloc(1,10,:,icarbon)',bm_alloc(1,10,:,icarbon),'frac_growthresp(10)',frac_growthresp(10)
467
468	!! 3.5 Total growth respiration
469	! Calculate total growth respiration of the plant as the sum of growth respiration of all plant parts
470	resp_growth(:,j) = zero
471
472	DO k = 1, nparts
473	resp_growth(:,j) = resp_growth(:,j) + resp_growth_part(:,k)
474	ENDDO
475
476	ENDDO ! # End Loop over # of PFTs
477	!WRITE(numout,*) 'zd nppcalc2 biomass(1,10,ileaf,icarbon)',biomass(1,10,ileaf,icarbon),'bm_pump(1)',bm_pump(1),'resp_maint_part(1,10,ileaf)',resp_maint_part(1,10,ileaf),'resp_maint(1,10)',resp_maint(1,10)
478	! WRITE(numout,*) 'frac_growthresp: ', frac_growthresp
479	! WRITE(numout,*) 'bm_alloc(1,12:14,:,icarbon) after update: ', bm_alloc(1,12:14,:,icarbon)
480
481	!! 4. Update the biomass with newly allocated biomass after respiration
482
483	! Save the old leaf biomass for later. "old" leaf mass is leaf mass after maintenance respiration in the case
484	! where maintenance respiration has required taking biomass from tissues in section 3.3
485	lm_old(:,:) = biomass(:,:,ileaf,icarbon)
486	biomass(:,:,:,:) = biomass(:,:,:,:) + bm_alloc(:,:,:,:)
487	!WRITE(numout,*) 'zd nppcalc3 biomass(1,10,ileaf,icarbon)',biomass(1,10,ileaf,icarbon)
488	!BEGINNVADD
489	npp_part (:,:,:)= npp_part (:,:,:) + bm_alloc(:,:,:,icarbon)
490	!ENDNVADD
491	! WRITE(numout,*) 'biomass(1,12:14,:,icarbon) after update: ', biomass(1,12:14,:,icarbon)
492	!! 5. Deal with negative biomasses
493
494	! Biomass can become negative in some rare cases, as the GPP can be negative. This corresponds to very
495	! situations that can be seen as the 'creation' of a seed ('virtual photosynthesis'). In this case, we set
496	! biomass to some small value min_stomate. For carbon budget to remain balanced, this creation of matter (carbon)
497	! is taken into account by decreasing the autotrophic respiration by the same amount that has been added to biomass
498	! for it to become positive. In this case, maintenance respiration can become negative in extreme cases (deserts)!!
499
500	DO k = 1, nparts ! Loop over # of plant parts
501
502	DO j = 2,nvm ! Loop over # of PFTs
503
504	WHERE ( biomass(:,j,k,icarbon) .LT. zero )
505
506	bm_create(:,j) = min_stomate - biomass(:,j,k,icarbon)
507
508	biomass(:,j,k,icarbon) = biomass(:,j,k,icarbon) + bm_create(:,j)
509
510	resp_maint(:,j) = resp_maint(:,j) - bm_create(:,j) / dt
511
512	ENDWHERE
513
514	ENDDO ! Loop over # of PFTs
515
516	ENDDO ! Loop over # plant parts
517	!WRITE(numout,*) 'zd nppcalc4 biomass(1,10,ileaf,icarbon)',biomass(1,10,ileaf,icarbon)
518
519	! WRITE(numout,*) 'biomass(1,12:14,:,icarbon) after negative removal: ', biomass(1,12:14,:,icarbon)
520	!! 6. Calculate NPP (See Eq 1 in header)
521
522	! Calculate the NPP @tex $(gC.m^{-2}dt^{-1})$ @endtex as the difference between GPP
523	! and autotrophic respiration (maintenance and growth respirations)
524	DO j = 2,nvm ! Loop over # PFTs
525	npp(:,j) = gpp(:,j) - resp_growth(:,j) - resp_maint(:,j)
526	ENDDO ! Loop over # PFTs
527
528
529	!! 7. Update leaf age
530
531	! Leaf age is needed for calculation of turnover and vmax in stomate_turnover.f90 and stomate_vmax.f90 routines.
532	! Leaf biomass is distributed according to its age into several "age classes" with age class=1 representing the
533	! youngest class, and consisting of the most newly allocated leaf biomass
534
535	!! 7.1 Update quantity and age of the leaf biomass in the youngest class
536	! The new amount of leaf biomass in the youngest age class (leaf_mass_young) is the sum of :
537	! - the leaf biomass that was already in the youngest age class (leaf_frac(:,j,1) * lm_old(:,j)) with the
538	! leaf age given in leaf_age(:,j,1)
539	! - and the new biomass allocated to leaves (bm_alloc(:,j,ileaf)) with a leaf age of zero.
540	DO j = 2,nvm
541	leaf_mass_young(:,j) = leaf_frac(:,j,1) * lm_old(:,j) + bm_alloc(:,j,ileaf,icarbon)
542	ENDDO
543
544	! The age of the updated youngest age class is the average of the ages of its 2 components: bm_alloc(leaf) of age
545	! '0', and leaf_frac*lm_old(=leaf_mass_young-bm_alloc) of age 'leaf_age(:,j,1)'
546	DO j = 2,nvm
547	WHERE ( ( bm_alloc(:,j,ileaf,icarbon) .GT. zero ) .AND. &
548	( leaf_mass_young(:,j) .GT. zero ) )
549
550	leaf_age(:,j,1) = MAX ( zero, &
551	& leaf_age(:,j,1) * &
552	& ( leaf_mass_young(:,j) - bm_alloc(:,j,ileaf,icarbon) ) / &
553	& leaf_mass_young(:,j) )
554
555	ENDWHERE
556	ENDDO
557
558	!! 7.2 Update leaf age
559	! Update fractions of leaf biomass in each age class (fraction in youngest class increases)
560
561	!! 7.2.1 Update age of youngest leaves
562	! For age class 1 (youngest class), because we have added biomass to the youngest class, we need to update
563	! the fraction of total leaf biomass that belongs to the youngest age class : updated mass in class divided
564	! by new total leaf mass
565	DO j = 2,nvm
566	WHERE ( biomass(:,j,ileaf,icarbon) .GT. min_stomate )
567
568	leaf_frac(:,j,1) = leaf_mass_young(:,j) / biomass(:,j,ileaf,icarbon)
569
570	ENDWHERE
571	ENDDO
572
573	!! 7.2.2 Update age of other age classes
574	! Because the total leaf biomass has changed, we need to update the fraction of leaves in each age class:
575	! mass in leaf age class (from previous fraction of leaves in this class and previous total leaf biomass)
576	! divided by new total mass
577	DO m = 2, nleafages ! Loop over # leaf age classes
578
579	DO j = 2,nvm ! Loop over # PFTs
580	WHERE ( biomass(:,j,ileaf,icarbon) .GT. min_stomate )
581
582	leaf_frac(:,j,m) = leaf_frac(:,j,m) * lm_old(:,j) / biomass(:,j,ileaf,icarbon)
583
584	ENDWHERE
585	ENDDO
586
587	ENDDO ! Loop over # leaf age classes
588	!gmjc varied sla for managed grassland
589	leaf_age_w = 0.0
590	DO j = 2,nvm
591	IF (is_grassland_manag(j)) THEN
592	sla_max_Nfert(:,j)=sla_max(j)
593	sla_min_Nfert(:,j)=sla_min(j)
594
595	ELSE
596	sla_max_Nfert(:,j)=sla_max(j)
597	sla_min_Nfert(:,j)=sla_min(j)
598	ENDIF
599
600	WHERE ( ( bm_alloc(:,j,ileaf,icarbon) .GT. 0.0 ) .AND. &
601	( leaf_mass_young(:,j) .GT. 0.0 ) )
602
603	sla_age1(:,j) = (sla_age1(:,j) * &
604	(leaf_mass_young(:,j)-bm_alloc(:,j,ileaf,icarbon)) + &
605	sla_max(j) * bm_alloc(:,j,ileaf,icarbon)) / leaf_mass_young(:,j)
606	sla_age2(:,j) = sla_max(j)*0.9
607	sla_age3(:,j) = sla_max(j)*0.85
608	sla_age4(:,j) = sla_max(j)*0.8
609
610	sla_calc(:,j) = sla_age1(:,j) * leaf_frac(:,j,1) + &
611	sla_age2(:,j) * leaf_frac(:,j,2) + &
612	sla_age3(:,j) * leaf_frac(:,j,3) + &
613	sla_age4(:,j) * leaf_frac(:,j,4)
614	ENDWHERE
615
616	leaf_age_w(:,j) = 0.0
617	DO m = 1, nleafages
618	leaf_age_w(:,j) = leaf_age_w(:,j)+ leaf_age(:,j,m)*leaf_frac(:,j,m)
619	END DO
620
621	! sla_calc can not be greater than sla_max or less than sla_min, and sla
622	! will be at maximum when age< 10
623	WHERE (sla_calc(:,j) .GT. sla_max(j))
624	sla_calc(:,j) = sla_max(j)
625	ELSE WHERE (leaf_age_w(:,j) .LT. 5.0)
626	sla_calc(:,j) = sla_max(j)
627	ELSE WHERE (sla_calc(:,j) .LT. sla_min(j))
628	sla_calc(:,j) = sla_min(j)
629	ENDWHERE
630	END DO
631	!end gmjc
632
633	!! 8. Update whole-plant age
634
635	!! 8.1 PFT age
636	! At every time step, increase age of the biomass that was already present at previous time step.
637	! Age is expressed in years, and the time step 'dt' in days so age increase is: dt divided by number
638	! of days in a year.
639	WHERE ( PFTpresent(:,:) )
640
641	age(:,:) = age(:,:) + dt/one_year
642
643	ELSEWHERE
644
645	age(:,:) = zero
646
647	ENDWHERE
648
649	!! 8.2 Age of grasses and crops
650	! For grasses and crops, biomass with age 0 has been added to the whole plant with age 'age'. New biomass is the sum of
651	! the current total biomass in all plant parts (bm_new), bm_new(:) = SUM( biomass(:,j,:), DIM=2 ). The biomass that has
652	! just been added is the sum of the allocatable biomass of all plant parts (bm_add), its age is zero. bm_add(:) =
653	! SUM( bm_alloc(:,j,:), DIM=2 ). Before allocation, the plant biomass is bm_new-bm_add, its age is "age(:,j)". The age of
654	! the new biomass is the average of the ages of previous and added biomass.
655	! For trees, age is treated in "establish" if vegetation is dynamic, and in turnover routines if it is static (in this
656	! case, only the age of the heartwood is accounted for).
657	DO j = 2,nvm
658
659	IF ( .NOT. is_tree(j) ) THEN
660
661	bm_new(:) = biomass(:,j,ileaf,icarbon) + biomass(:,j,isapabove,icarbon) + &
662	biomass(:,j,iroot,icarbon) + biomass(:,j,ifruit,icarbon)
663	bm_add(:) = bm_alloc(:,j,ileaf,icarbon) + bm_alloc(:,j,isapabove,icarbon) + &
664	bm_alloc(:,j,iroot,icarbon) + bm_alloc(:,j,ifruit,icarbon)
665
666	WHERE ( ( bm_new(:) .GT. zero ) .AND. ( bm_add(:) .GT. zero ) )
667	age(:,j) = age(:,j) * ( bm_new(:) - bm_add(:) ) / bm_new(:)
668	ENDWHERE
669
670	ENDIF
671
672	ENDDO
673
674	!! 9. Write history files
675
676	CALL xios_orchidee_send_field("BM_ALLOC_LEAF",bm_alloc(:,:,ileaf,icarbon))
677	CALL xios_orchidee_send_field("BM_ALLOC_SAP_AB",bm_alloc(:,:,isapabove,icarbon))
678	CALL xios_orchidee_send_field("BM_ALLOC_SAP_BE",bm_alloc(:,:,isapbelow,icarbon))
679	CALL xios_orchidee_send_field("BM_ALLOC_ROOT",bm_alloc(:,:,iroot,icarbon))
680	CALL xios_orchidee_send_field("BM_ALLOC_FRUIT",bm_alloc(:,:,ifruit,icarbon))
681	CALL xios_orchidee_send_field("BM_ALLOC_RES",bm_alloc(:,:,icarbres,icarbon))
682	!gmjc
683	npp_above(:,:) = npp_part(:,:,ileaf)+npp_part(:,:,isapabove)+&
684	npp_part(:,:,ifruit)+npp_part(:,:,icarbres)/2.
685	npp_below(:,:) = npp_part(:,:,iroot)+npp_part(:,:,isapbelow)+&
686	npp_part(:,:,icarbres)/2.
687	CALL xios_orchidee_send_field("SLA_CALC",sla_calc)
688	CALL xios_orchidee_send_field("NPP_ABOVE",npp_above)
689	CALL xios_orchidee_send_field("NPP_BELOW",npp_below)
690	!end gmjc
691
692	! Save in history file the variables describing the biomass allocated to the plant parts
693	CALL histwrite_p (hist_id_stomate, 'SLA_CROP',itime, &
694	ssla, npts*nvm, horipft_index)
695	CALL histwrite_p (hist_id_stomate, 'BM_ALLOC_LEAF', itime, &
696	bm_alloc(:,:,ileaf,icarbon), npts*nvm, horipft_index)
697	CALL histwrite_p (hist_id_stomate, 'BM_ALLOC_SAP_AB', itime, &
698	bm_alloc(:,:,isapabove,icarbon), npts*nvm, horipft_index)
699	CALL histwrite_p (hist_id_stomate, 'BM_ALLOC_SAP_BE', itime, &
700	bm_alloc(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
701	CALL histwrite_p (hist_id_stomate, 'BM_ALLOC_ROOT', itime, &
702	bm_alloc(:,:,iroot,icarbon), npts*nvm, horipft_index)
703	CALL histwrite_p (hist_id_stomate, 'BM_ALLOC_FRUIT', itime, &
704	bm_alloc(:,:,ifruit,icarbon), npts*nvm, horipft_index)
705	CALL histwrite_p (hist_id_stomate, 'BM_ALLOC_RES', itime, &
706	bm_alloc(:,:,icarbres,icarbon), npts*nvm, horipft_index)
707	!gmjc
708	CALL histwrite_p(hist_id_stomate, 'SLA_CALC', itime, &
709	sla_calc(:,:), npts*nvm, horipft_index)
710	!end gmjc
711	!BEGINNVADD
712	CALL histwrite_p(hist_id_stomate, 'NPP_ABOVE', itime, &
713	npp_above, npts*nvm, horipft_index)
714	CALL histwrite_p(hist_id_stomate, 'NPP_BELOW', itime, &
715	npp_below, npts*nvm, horipft_index)
716	!ENDNVADD
717	IF (printlev>=4) WRITE(numout,*) 'Leaving npp'
718
719	END SUBROUTINE npp_calc
720
721	END MODULE stomate_npp

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd_mict_peat_ch4/src_stomate/stomate_npp.f90 @ 7346

Download in other formats: