Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_stomate/stomate_litter.f90 @ 5491

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

Line
1	! =================================================================================================================================
2	! MODULE : stomate_litter
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 Update litter and lignine content after litter fall and
10	!! calculating litter decomposition.
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! REFERENCE(S) : None
17	!!
18	!! SVN :
19	!! $HeadURL$
20	!! $Date$
21	!! $Revision$
22	!! \n
23	!_ ================================================================================================================================
24
25	MODULE stomate_litter
26
27	! modules used:
28
29	USE ioipsl_para
30	USE stomate_data
31	USE constantes
32	USE constantes_soil
33	USE pft_parameters
34
35	IMPLICIT NONE
36
37	! private & public routines
38
39	PRIVATE
40	PUBLIC littercalc,littercalc_clear, deadleaf
41
42	LOGICAL, SAVE :: firstcall_litter = .TRUE. !! first call
43	!$OMP THREADPRIVATE(firstcall_litter)
44
45	CONTAINS
46
47	!! ================================================================================================================================
48	!! SUBROUTINE : littercalc_calc
49	!!
50	!!\BRIEF Set the flag ::firstcall_litter to .TRUE. and as such activate section
51	!! 1.1 of the subroutine littercalc (see below).
52	!!
53	!! DESCRIPTION : None
54	!!
55	!! RECENT CHANGE(S) : None
56	!!
57	!! MAIN OUTPUT VARIABLE(S) : None
58	!!
59	!! REFERENCE(S) : None
60	!!
61	!! FLOWCHART : None
62	!! \n
63	!_ ================================================================================================================================
64
65	SUBROUTINE littercalc_clear
66	firstcall_litter =.TRUE.
67	END SUBROUTINE littercalc_clear
68
69
70	!! ================================================================================================================================
71	!! SUBROUTINE : littercalc
72	!!
73	!!\BRIEF Calculation of the litter decomposition and therefore of the
74	!! heterotrophic respiration from litter following Parton et al. (1987).
75	!!
76	!! DESCRIPTION : The littercal routine splits the litter in 4 pools:
77	!! aboveground metaboblic, aboveground structural, belowground metabolic and
78	!! belowground structural. the fraction (F) of plant material going to metabolic
79	!! and structural is defined following Parton et al. (1987)
80	!! \latexonly
81	!! \input{littercalc1.tex}
82	!! \endlatexonly
83	!! \n
84	!! where L is the lignin content of the plant carbon pools considered and CN
85	!! its CN ratio. L and CN are fixed parameters for each plant carbon pools,
86	!! therefore it is the ratio between each plant carbon pool within a PFT, which
87	!! controlled the part of the total litter, that will be considered as
88	!! recalcitrant (i.e. structural litter) or labile (i.e. metabolic litter).\n
89	!!
90	!! The routine calculates the fraction of aboveground litter which is metabolic
91	!! or structural (the litterpart variable) which is then used in lpj_fire.f90.\n
92	!!
93	!! In the section 2, the routine calculate the new plant material entering the
94	!! litter pools by phenological death of plants organs (corresponding to the
95	!! variable turnover) and by fire, herbivory and others non phenological causes
96	!! (variable bm_to_litter). This calculation is first done for each PFT and then
97	!! the values calculated for each PFT are added up. Following the same approach
98	!! the lignin content of the total structural litter is calculated and will be
99	!! then used as a factor control of the decomposition of the structural litter
100	!! (lignin_struc) in the section 5.1.2. A test is performed to avoid that we add
101	!! more lignin than structural litter. Finally, the variable litterpart is
102	!! updated.\n
103	!!
104	!! In the section 3 and 4 the temperature and the moisture controlling the
105	!! decomposition are calculated for above and belowground. For aboveground
106	!! litter, air temperature and litter moisture are calculated in sechiba and used
107	!! directly. For belowground, soil temperature and moisture are also calculated
108	!! in sechiba but are modulated as a function of the soil depth. The modulation
109	!! is a multiplying factor exponentially distributed between 0 (in depth) and 1
110	!! in surface.\n
111	!!
112	!! Then, in the section 5, the routine calculates the structural litter decomposition
113	!! (C) following first order kinetics following Parton et al. (1987).
114	!! \latexonly
115	!! \input{littercalc2.tex}
116	!! \endlatexonly
117	!! \n
118	!! with k the decomposition rate of the structural litter.
119	!! k corresponds to
120	!! \latexonly
121	!! \input{littercalc3.tex}
122	!! \endlatexonly
123	!! \n
124	!! with littertau the turnover rate, T a function of the temperature and M a function of
125	!! the moisture described below.\n
126	!!
127	!! Then, the fraction of dead leaves (DL) composed by aboveground structural litter is
128	!! calculated as following
129	!! \latexonly
130	!! \input{littercalc4.tex}
131	!! \endlatexonly
132	!! \n
133	!! with k the decomposition rate of the structural litter previously
134	!! described.\n
135	!!
136	!! In the section 5.1, the fraction of decomposed structural litter
137	!! incorporated to the soil (Input) and its associated heterotrophic respiration are
138	!! calculated. For structural litter, the C decomposed could go in the active
139	!! soil carbon pool or in the slow carbon, as described in
140	!! stomate_soilcarbon.f90.\n
141	!! \latexonly
142	!! \input{littercalc5.tex}
143	!! \endlatexonly
144	!! \n
145	!! with f a parameter describing the fraction of structural litter incorporated
146	!! into the considered soil carbon pool, C the amount of litter decomposed and L
147	!! the amount of lignin in the litter. The litter decomposed which is not
148	!! incorporated into the soil is respired.\n
149	!!
150	!! In the section 5.2, the fraction of decomposed metabolic litter
151	!! incorporated to the soil and its associated heterotrophic respiration are
152	!! calculated with the same approaches presented for 5.1 but no control factor
153	!! depending on the lignin content are used.\n
154	!!
155	!! In the section 6 the dead leaf cover is calculated through a call to the
156	!! deadleaf subroutine presented below.\n
157	!!
158	!! In the section 7, if the flag SPINUP_ANALYTIC is set to true, we fill MatrixA
159	!! and VectorB following Lardy(2011).
160	!!
161	!! MAIN OUTPUT VARIABLES: ::deadleaf_cover, ::resp_hetero_litter, ::soilcarbon_input,
162	!! ::control_temp, ::control_moist
163	!!
164	!! REFERENCES:
165	!! - Parton, WJ, Schimel, DS, Cole, CV, and Ojima, DS. 1987. Analysis
166	!! of factors controlling soil organic matter levels in Great Plains
167	!! grasslands. Soil Science Society of America journal (USA)
168	!! (51):1173-1179.
169	!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
170	!! Environmental Modelling & Software (2011), doi:10.1016\|j.envsoft.2011.05.016
171	!!
172	!! FLOWCHART :
173	!! \latexonly
174	!! \includegraphics(scale=0.5){littercalcflow.jpg}
175	!! \endlatexonly
176	!! \n
177	!_ ================================================================================================================================
178
179	SUBROUTINE littercalc (npts, &
180	turnover, bm_to_litter, &
181	veget_max, tsurf, tsoil, soilhum, litterhum, &
182	litterpart, litter, litter_avail, litter_not_avail, litter_avail_frac, &
183	!spitfire
184	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
185	!endspit
186	dead_leaves, lignin_struc, &
187	deadleaf_cover, resp_hetero_litter, &
188	soilcarbon_input, control_temp, control_moist, &
189	MatrixA, VectorB, &!)
190	!!!!!crops
191	tsoil_cm, soilhum_cm, &
192	!!!!!xuhui
193	!gmjc
194	sla_calc,do_slow, & !)
195	!end gmjc
196	!!!qcj++ peatland
197	shumdiag_peat,mc_peat_above)
198
199	!! 0. Variable and parameter declaration
200
201	!! 0.1 Input variables
202
203	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixels
204	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: turnover !! Turnover rates of plant biomass
205	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
206	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: bm_to_litter !! Conversion of biomass to litter
207	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
208	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
209	!! defined in the input vegetation map
210	!! @tex $(m^2 m^{-2})$ @endtex
211	REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf !! Temperature (K) at the surface
212	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: tsoil !! Soil temperature (K)
213	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: soilhum !! Daily soil humidity of each soil layer
214	!! (unitless)
215	REAL(r_std), DIMENSION(npts), INTENT(in) :: litterhum !! Daily litter humidity (unitless)
216
217	!! 0.2 Output variables
218
219	REAL(r_std), DIMENSION(npts), INTENT(out) :: deadleaf_cover !! Fraction of soil covered by dead leaves
220	!! over all PFTs (0-1, unitless)
221	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_hetero_litter !! Litter heterotrophic respiration. The unit
222	!! is given by m^2 of ground.
223	!! @tex $(gC dt_sechiba one\_day^{-1}) m^{-2})$ @endtex
224	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(out) :: soilcarbon_input !! Quantity of carbon going into carbon pools
225	!! from litter decomposition. The unit is
226	!! given by m^2 of ground
227	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
228	REAL(r_std), DIMENSION(npts,nlevs), INTENT(out) :: control_temp !! Temperature control of heterotrophic
229	!! respiration, above and below (0-1,
230	!! unitless)
231	REAL(r_std), DIMENSION(npts,nlevs), INTENT(out) :: control_moist !! Moisture control of heterotrophic
232	!! respiration (0.25-1, unitless)
233	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(out) :: MatrixA !! Matrix containing the fluxes between the
234	!! carbon pools per sechiba time step
235	!! @tex $(gC.m^2.day^{-1})$ @endtex
236	REAL(r_std), DIMENSION(npts,nvm,nbpools), INTENT(out) :: VectorB !! Vector containing the litter increase per
237	!! sechiba time step
238	!! @tex $(gC m^{-2})$ @endtex
239	!!!qcj++ peatland
240	REAL(r_std),DIMENSION (npts,nslm),INTENT(in) :: shumdiag_peat
241	REAL(r_std), DIMENSION(npts,nlevs) :: control_moist_peat
242	REAL(r_std), DIMENSION(npts,nlevs) :: control_temp_peat
243	REAL(r_std), DIMENSION(npts),INTENT(in) :: mc_peat_above
244	REAL(r_std), DIMENSION(npts) :: mc_peat_below
245
246	!!!!! crops
247	! soil temperature at the resolution of 1 cm, the second dimension is 3
248	! according to the STICS code. That represent the 3 layers around sowing
249	! depth
250	REAL(r_std), DIMENSION(npts, nvm, 3), INTENT(inout) :: tsoil_cm
251	! soil relative humidity at the resolution of 1 cm, the second dimension is
252	! 3
253	! according to the STICS code. That represent the 3 layers around sowing
254	! depth
255	REAL(r_std), DIMENSION(npts, nvm, 3), INTENT(inout) :: soilhum_cm
256	!!!!! end crops, xuhui
257
258	!! 0.3 Modified variables
259
260	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: litterpart !! Fraction of litter above the ground
261	!! belonging to the different PFTs (0-1,
262	!! unitless)
263	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: litter !! Metabolic and structural litter,above and
264	!! below ground. The unit is given by m^2 of
265	!! ground @tex $(gC m^{-2})$ @endtex
266	REAL(r_std), DIMENSION(npts,nlitt,nvm), INTENT(inout) :: litter_avail
267	REAL(r_std), DIMENSION(npts,nlitt,nvm), INTENT(in) :: litter_not_avail
268	REAL(r_std), DIMENSION(npts,nlitt,nvm), INTENT(out) :: litter_avail_frac
269
270	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: dead_leaves !! Dead leaves per ground unit area, per PFT,
271	!! metabolic and structural in
272	!! @tex $(gC m^{-2})$ @endtex
273	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout) :: lignin_struc !! Ratio Lignin content in structural litter,
274	!! above and below ground, (0-1, unitless)
275	!gmjc
276	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: sla_calc
277	LOGICAL,INTENT(in) :: do_slow
278	!end gmjc
279	!! 0.4 Local variables
280
281	REAL(r_std) :: dt !! Time step of the simulations for stomate
282	!! @tex $(dt_sechiba one\_day^{-1})$ @endtex
283	REAL(r_std), SAVE, DIMENSION(nparts,nlitt) :: litterfrac !! The fraction of leaves, wood, etc. that
284	!! goes into metabolic and structural
285	!! litterpools (0-1, unitless)
286	!$OMP THREADPRIVATE(litterfrac)
287	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: z_soil !! Soil levels (m)
288	!$OMP THREADPRIVATE(z_soil)
289	REAL(r_std), DIMENSION(npts) :: rpc !! Integration constant for vertical root
290	!! profiles (unitless)
291	REAL(r_std), SAVE, DIMENSION(nlitt) :: litter_tau !! Turnover time in litter pools (days)
292	!$OMP THREADPRIVATE(litter_tau)
293	REAL(r_std), SAVE, DIMENSION(nlitt,ncarb,nlevs) :: frac_soil !! Fraction of litter that goes into soil
294	!! (litter -> carbon, above and below). The
295	!! remaining part goes to the atmosphere
296	!$OMP THREADPRIVATE(frac_soil)
297	REAL(r_std), DIMENSION(npts) :: tsoil_decomp !! Temperature used for decompostition in
298	!! soil (K)
299	REAL(r_std), DIMENSION(npts) :: soilhum_decomp !! Humidity used for decompostition in soil
300	!! (unitless)
301	REAL(r_std), DIMENSION(npts) :: fd !! Fraction of structural or metabolic litter
302	!! decomposed (unitless)
303	REAL(r_std), DIMENSION(npts,nelements) :: qd !! Quantity of structural or metabolic litter
304	!! decomposed @tex $(gC m^{-2})$ @endtex
305	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: old_struc !! Old structural litter, above and below
306	!! @tex $(gC m^{-2})$ @endtex
307	REAL(r_std), DIMENSION(npts,nvm,nlitt,nlevs,nelements) :: litter_inc_PFT !! Increase of litter, per PFT, metabolic and
308	!! structural, above and below ground. The
309	!! unit is given by m^2 of ground.
310	!! @tex $(gC m^{-2})$ @endtex
311	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements) :: litter_inc !! Increase of metabolic and structural
312	!! litter, above and below ground. The unit
313	!! is given by m^2 of ground.
314	!! @tex $(gC m^{-2})$ @endtex
315	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: lignin_struc_inc !! Lignin increase in structural litter,
316	!! above and below ground. The unit is given
317	!! by m^2 of ground.
318	!! @tex $(gC m^{-2})$ @endtex
319	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements) :: litter_pft !! Metabolic and structural litter above the
320	!! ground per PFT
321	REAL(r_std), DIMENSION(npts) :: zdiff_min !! Intermediate field for looking for minimum
322	!! of what? this is not used in the code.
323	!! [??CHECK] could we delete it?
324	CHARACTER(LEN=10), DIMENSION(nlitt) :: litter_str !! Messages to write output information about
325	!! the litter
326	CHARACTER(LEN=22), DIMENSION(nparts) :: part_str !! Messages to write output information about
327	!! the plant
328	CHARACTER(LEN=7), DIMENSION(ncarb) :: carbon_str !! Messages to write output information about
329	!! the soil carbon
330	CHARACTER(LEN=5), DIMENSION(nlevs) :: level_str !! Messages to write output information about
331	!! the level (aboveground or belowground litter)
332	INTEGER(i_std) :: i,j,k,l,m !! Indices (unitless)
333	INTEGER(i_std), SAVE :: frozen_respiration_func = 1
334	LOGICAL, SAVE :: use_snowinsul_litter_resp = .FALSE.
335
336	!spitfire
337	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_1hr
338	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_10hr
339	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_100hr
340	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_1000hr
341	!total fuel load summing together all fuel classes.
342	REAl(r_std), DIMENSION(npts,nvm,nlitt,nelements) :: fuel_nlitt_total_pft
343	! difference between litter and fuel_nlitt_total_pft
344	REAL(r_std), DIMENSION(npts) :: diff_frac
345	! quantity of structural or metabolic litter decomposed (gC/m**2) for decomposing fuel classes per PFT
346	! chaocomm: this is now only a temporary variable and will be deleted later.
347	REAL(r_std), DIMENSION(npts,nelements) :: qd_fuel
348	!endspit
349	!!!!! crops
350	! local variables
351	! three layer around the sowing depth, in cm
352	REAL(r_std), DIMENSION(nvm, 3) :: ldep
353
354	! the nth of half layer,[1, 2, .....2*nbdl]
355	INTEGER(i_std) :: nhl
356	! the upper limit humidity for ith half layer of 22 half layers, 2*nbdl
357	INTEGER(i_std), DIMENSION(npts, 2*nslm) :: ulhum
358	! the ith layer for the three layers around sowing depth
359	INTEGER(i_std) :: il
360	! uldep is the upper limit of ith half layer of 22 half layers, 2*nbdl
361	!REAL(r_std), DIMENSION(2*nbdl) :: uldep
362	REAL(r_std) :: uldep
363	! bldep is the lower limit of ith half layer of 22 half layers
364	REAL(r_std) :: bldep
365
366	! the upper limit soil temperature for ith half layer of 22 half layers, 2*nbdl
367	REAL(r_std), DIMENSION(npts, 2*nslm) :: ultsoil
368	! identifier for number of PFTs
369	INTEGER(i_std) :: ipft
370	!!!!! end crops
371
372	INTEGER(i_std) :: ier !! Error handling
373	!_ ================================================================================================================================
374
375	IF (printlev>=3) WRITE(numout,*) 'Entering littercalc'
376
377	!! 1. Initialisations of the different fields during the first call of the routine
378	dt = dt_sechiba/one_day
379	diff_frac(:) = 0
380
381	IF ( firstcall_litter ) THEN
382
383	!! 1.1.3 litter fractions:
384	!! what fraction of leaves, wood, etc. goes into metabolic and structural litterpools
385	DO k = 1, nparts
386
387	litterfrac(k,imetabolic) = metabolic_ref_frac - metabolic_LN_ratio * LC(k) * CN(k)
388	litterfrac(k,istructural) = un - litterfrac(k,imetabolic)
389
390	ENDDO
391
392	!! 1.1.4 residence times in litter pools (days)
393	litter_tau(imetabolic) = tau_metabolic * one_year ! .5 years
394	litter_tau(istructural) = tau_struct * one_year ! 3 years
395
396	!! 1.1.5 decomposition flux fraction that goes into soil
397	! (litter -> carbon, above and below)
398	! 1-frac_soil goes into atmosphere
399	frac_soil(:,:,:) = zero
400
401	! structural litter: lignin fraction goes into slow pool + respiration,
402	! rest into active pool + respiration
403	frac_soil(istructural,iactive,iabove) = frac_soil_struct_aa
404	frac_soil(istructural,iactive,ibelow) = frac_soil_struct_ab
405	frac_soil(istructural,islow,iabove) = frac_soil_struct_sa
406	frac_soil(istructural,islow,ibelow) = frac_soil_struct_sb
407
408	! metabolic litter: all goes into active pool + respiration.
409	! Nothing into slow or passive pool.
410	frac_soil(imetabolic,iactive,iabove) = frac_soil_metab_aa
411	frac_soil(imetabolic,iactive,ibelow) = frac_soil_metab_ab
412
413
414	!! 1.2 soil levels
415	ALLOCATE(z_soil(0:nslm), stat=ier)
416	IF ( ier /= 0 ) CALL ipslerr_p(3,'littercalc','Pb in allocate of z_soil','','')
417
418	z_soil(0) = zero
419	z_soil(1:nslm) = diaglev(1:nslm)
420
421
422	!! 1.3 messages
423	litter_str(imetabolic) = 'metabolic'
424	litter_str(istructural) = 'structural'
425
426	carbon_str(iactive) = 'active'
427	carbon_str(islow) = 'slow'
428	carbon_str(ipassive) = 'passive'
429
430	level_str(iabove) = 'above'
431	level_str(ibelow) = 'below'
432
433	part_str(ileaf) = 'leaves'
434	part_str(isapabove) = 'sap above ground'
435	part_str(isapbelow) = 'sap below ground'
436	part_str(iheartabove) = 'heartwood above ground'
437	part_str(iheartbelow) = 'heartwood below ground'
438	part_str(iroot) = 'roots'
439	part_str(ifruit) = 'fruits'
440	part_str(icarbres) = 'carbohydrate reserve'
441
442	WRITE(numout,*) 'litter:'
443
444	WRITE(numout,*) ' > C/N ratios: '
445	DO k = 1, nparts
446	WRITE(numout,*) ' ', part_str(k), ': ',CN(k)
447	ENDDO
448
449	WRITE(numout,*) ' > Lignine/C ratios: '
450	DO k = 1, nparts
451	WRITE(numout,*) ' ', part_str(k), ': ',LC(k)
452	ENDDO
453
454	WRITE(numout,*) ' > fraction of compartment that goes into litter: '
455	DO k = 1, nparts
456	DO m = 1, nlitt
457	WRITE(numout,*) ' ', part_str(k), '-> ',litter_str(m), ':',litterfrac(k,m)
458	ENDDO
459	ENDDO
460
461	WRITE(numout,*) ' > scaling depth for decomposition (m): ',z_decomp
462
463	WRITE(numout,*) ' > minimal carbon residence time in litter pools (d):'
464	DO m = 1, nlitt
465	WRITE(numout,*) ' ',litter_str(m),':',litter_tau(m)
466	ENDDO
467
468	WRITE(numout,*) ' > litter decomposition flux fraction that really goes '
469	WRITE(numout,*) ' into carbon pools (rest into the atmosphere):'
470	DO m = 1, nlitt
471	DO l = 1, nlevs
472	DO k = 1, ncarb
473	WRITE(numout,*) ' ',litter_str(m),' ',level_str(l),' -> ',&
474	carbon_str(k),':', frac_soil(m,k,l)
475	ENDDO
476	ENDDO
477	ENDDO
478
479	CALL getin_p('frozen_respiration_func',frozen_respiration_func)
480	WRITE(numout, *)' frozen litter respiration function: ', frozen_respiration_func
481
482	CALL getin_p('use_snowinsul_litter_resp',use_snowinsul_litter_resp)
483	WRITE(numout, *)' use_snowinsul_litter_resp ( i.e. set surface litter temp to below snow layer:) ', use_snowinsul_litter_resp
484
485	firstcall_litter = .FALSE.
486
487	ENDIF
488
489
490	!! 1.3 litter above the ground per PFT.
491	DO j = 1, nvm ! Loop over # PFTs
492
493	DO k = 1, nlitt !Loop over litter pool
494	litter_pft(:,j,k,icarbon) = litterpart(:,j,k) * litter(:,k,j,iabove,icarbon)
495	ENDDO
496
497	ENDDO
498
499
500	!! 1.4 set output to zero
501	deadleaf_cover(:) = zero
502	resp_hetero_litter(:,:) = zero
503	soilcarbon_input(:,:,:) = zero
504
505
506	!! 2. Add biomass to different litterpools (per m^2 of ground)
507
508	!! 2.1 first, save old structural litter (needed for lignin fractions).
509	! above/below
510	DO l = 1, nlevs !Loop over litter levels (above and below ground)
511	DO m = 1,nvm !Loop over PFTs
512
513	old_struc(:,m,l) = litter(:,istructural,m,l,icarbon)
514
515	ENDDO
516	ENDDO
517
518
519	!! 2.2 update litter, dead leaves, and lignin content in structural litter
520	litter_inc(:,:,:,:,:) = zero
521	lignin_struc_inc(:,:,:) = zero
522
523	DO j = 1,nvm !Loop over PFTs
524
525	!! 2.2.1 litter
526	DO k = 1, nlitt !Loop over litter pools (metabolic and structural)
527
528	!! 2.2.2 calculate litter increase (per m^2 of ground).
529	! Only a given fracion of fruit turnover is directly coverted into litter.
530	! Litter increase for each PFT, structural and metabolic, above/below
531	litter_inc_PFT(:,j,k,iabove,icarbon) = &
532	litterfrac(ileaf,k) * bm_to_litter(:,j,ileaf,icarbon) + &
533	litterfrac(isapabove,k) * bm_to_litter(:,j,isapabove,icarbon) + &
534	litterfrac(iheartabove,k) * bm_to_litter(:,j,iheartabove,icarbon) + &
535	litterfrac(ifruit,k) * bm_to_litter(:,j,ifruit,icarbon) + &
536	litterfrac(icarbres,k) * bm_to_litter(:,j,icarbres,icarbon) + &
537	litterfrac(ileaf,k) * turnover(:,j,ileaf,icarbon) + &
538	litterfrac(isapabove,k) * turnover(:,j,isapabove,icarbon) + &
539	litterfrac(iheartabove,k) * turnover(:,j,iheartabove,icarbon) + &
540	litterfrac(ifruit,k) * turnover(:,j,ifruit,icarbon) + &
541	litterfrac(icarbres,k) * turnover(:,j,icarbres,icarbon)
542
543	litter_inc_PFT(:,j,k,ibelow,icarbon) = &
544	litterfrac(isapbelow,k) * bm_to_litter(:,j,isapbelow,icarbon) + &
545	litterfrac(iheartbelow,k) * bm_to_litter(:,j,iheartbelow,icarbon) + &
546	litterfrac(iroot,k) * bm_to_litter(:,j,iroot,icarbon) + &
547	litterfrac(isapbelow,k) * turnover(:,j,isapbelow,icarbon) + &
548	litterfrac(iheartbelow,k) * turnover(:,j,iheartbelow,icarbon) + &
549	litterfrac(iroot,k) * turnover(:,j,iroot,icarbon)
550
551	! litter increase, met/struct, above/below
552	litter_inc(:,k,j,iabove,icarbon) = litter_inc(:,k,j,iabove,icarbon) + litter_inc_PFT(:,j,k,iabove,icarbon)
553	litter_inc(:,k,j,ibelow,icarbon) = litter_inc(:,k,j,ibelow,icarbon) + litter_inc_PFT(:,j,k,ibelow,icarbon)
554
555	!! 2.2.3 dead leaves, for soil cover.
556	dead_leaves(:,j,k) = &
557	dead_leaves(:,j,k) + &
558	litterfrac(ileaf,k) * ( bm_to_litter(:,j,ileaf,icarbon) + turnover(:,j,ileaf,icarbon) )
559
560	!! 2.2.4 lignin increase in structural litter
561	IF ( k .EQ. istructural ) THEN
562
563	lignin_struc_inc(:,j,iabove) = &
564	lignin_struc_inc(:,j,iabove) + &
565	LC(ileaf) * bm_to_litter(:,j,ileaf,icarbon) + &
566	LC(isapabove) * bm_to_litter(:,j,isapabove,icarbon) + &
567	LC(iheartabove) * bm_to_litter(:,j,iheartabove,icarbon) + &
568	LC(ifruit) * bm_to_litter(:,j,ifruit,icarbon) + &
569	LC(icarbres) * bm_to_litter(:,j,icarbres,icarbon) + &
570	LC(ileaf) * turnover(:,j,ileaf,icarbon) + &
571	LC(isapabove) * turnover(:,j,isapabove,icarbon) + &
572	LC(iheartabove) * turnover(:,j,iheartabove,icarbon) + &
573	LC(ifruit) * turnover(:,j,ifruit,icarbon) + &
574	LC(icarbres) * turnover(:,j,icarbres,icarbon)
575
576	lignin_struc_inc(:,j,ibelow) = &
577	lignin_struc_inc(:,j,ibelow) + &
578	LC(isapbelow) * bm_to_litter(:,j,isapbelow,icarbon) + &
579	LC(iheartbelow) * bm_to_litter(:,j,iheartbelow,icarbon) + &
580	LC(iroot) * bm_to_litter(:,j,iroot,icarbon) + &
581	LC(isapbelow)*turnover(:,j,isapbelow,icarbon) + &
582	LC(iheartbelow)*turnover(:,j,iheartbelow,icarbon) + &
583	LC(iroot)*turnover(:,j,iroot,icarbon)
584
585	ENDIF
586
587	!spitfire
588	!calculate fuel amount
589	fuel_1hr(:,j,k,icarbon) = fuel_1hr(:,j,k,icarbon) + &
590	1.0 * (litterfrac(ileaf,k) * bm_to_litter(:,j,ileaf,icarbon) + &
591	litterfrac(ileaf,k) * turnover(:,j,ileaf,icarbon) + &
592	litterfrac(ifruit,k) * bm_to_litter(:,j,ifruit,icarbon) + &
593	litterfrac(ifruit,k) * turnover(:,j,ifruit,icarbon)) + &
594	0.0450 * (litterfrac(isapabove,k) * bm_to_litter(:,j,isapabove,icarbon) + &
595	litterfrac(isapabove,k) * turnover(:,j,isapabove,icarbon) + &
596	litterfrac(iheartabove,k) * bm_to_litter(:,j,iheartabove,icarbon) + &
597	litterfrac(iheartabove,k) * turnover(:,j,iheartabove,icarbon) + &
598	litterfrac(icarbres,k) * bm_to_litter(:,j,icarbres,icarbon) + &
599	litterfrac(icarbres,k) * turnover(:,j,icarbres,icarbon) )
600
601
602	fuel_10hr(:,j,k,icarbon) = fuel_10hr(:,j,k,icarbon) + &
603	0.0750 * (litterfrac(isapabove,k) * bm_to_litter(:,j,isapabove,icarbon) + &
604	litterfrac(isapabove,k) * turnover(:,j,isapabove,icarbon) + &
605	litterfrac(iheartabove,k) * bm_to_litter(:,j,iheartabove,icarbon) + &
606	litterfrac(iheartabove,k) * turnover(:,j,iheartabove,icarbon) + &
607	litterfrac(icarbres,k) * bm_to_litter(:,j,icarbres,icarbon) + &
608	litterfrac(icarbres,k) * turnover(:,j,icarbres,icarbon))
609
610	fuel_100hr(:,j,k,icarbon) = fuel_100hr(:,j,k,icarbon) + &
611	0.2100 * (litterfrac(isapabove,k) * bm_to_litter(:,j,isapabove,icarbon) + &
612	litterfrac(isapabove,k) * turnover(:,j,isapabove,icarbon) + &
613	litterfrac(iheartabove,k) * bm_to_litter(:,j,iheartabove,icarbon) + &
614	litterfrac(iheartabove,k) * turnover(:,j,iheartabove,icarbon) + &
615	litterfrac(icarbres,k) * bm_to_litter(:,j,icarbres,icarbon) + &
616	litterfrac(icarbres,k) * turnover(:,j,icarbres,icarbon))
617
618	fuel_1000hr(:,j,k,icarbon) = fuel_1000hr(:,j,k,icarbon) + &
619	0.6700 * (litterfrac(isapabove,k) * bm_to_litter(:,j,isapabove,icarbon) + &
620	litterfrac(isapabove,k) * turnover(:,j,isapabove,icarbon) + &
621	litterfrac(iheartabove,k) * bm_to_litter(:,j,iheartabove,icarbon) + &
622	litterfrac(iheartabove,k) * turnover(:,j,iheartabove,icarbon) + &
623	litterfrac(icarbres,k) * bm_to_litter(:,j,icarbres,icarbon) + &
624	litterfrac(icarbres,k) * turnover(:,j,icarbres,icarbon))
625
626	fuel_nlitt_total_pft(:,j,k,icarbon) = fuel_1hr(:,j,k,icarbon) + fuel_10hr(:,j,k,icarbon) + &
627	fuel_100hr(:,j,k,icarbon) + fuel_1000hr(:,j,k,icarbon)
628
629	!endspit
630	ENDDO
631	ENDDO
632
633	!! 2.2.5 add new litter (struct/met, above/below)
634	litter(:,:,:,:,:) = litter(:,:,:,:,:) + litter_inc(:,:,:,:,:)
635	!gmjc for grazing litter
636	IF (do_slow) THEN
637	DO j=2,nvm
638	!! grazed grassland or natural grassland with wild animal
639	IF ((is_grassland_manag(j) .AND. is_grassland_grazed(j)) .OR. &
640	((.NOT. is_tree(j)) .AND. natural(j) .AND. &
641	& (.NOT. is_grassland_cut(j)) .AND. (.NOT.is_grassland_grazed(j)))) THEN
642	WHERE (litter(:,:,j,iabove,icarbon) .GE. litter_not_avail(:,:,j) &
643	.AND. litter(:,:,j,iabove,icarbon) .GT. 0.0)
644	litter_avail_frac(:,:,j) = (litter(:,:,j,iabove,icarbon) - litter_not_avail(:,:,j)) &
645	& / litter(:,:,j,iabove,icarbon)
646	! if litter not available equal to or larger than litter
647	ELSEWHERE
648	litter_avail_frac(:,:,j) = 0.0
649	ENDWHERE
650	IF (ANY(litter_avail_frac .LT. 0.0)) THEN
651	WRITE (numout,*) 'frac error',litter(:,:,j,iabove,icarbon),litter_not_avail(:,:,j),litter_avail_frac(:,:,j)
652	STOP 'error fraction litter available for grazing < 0'
653	ENDIF
654	ELSE
655	litter_avail_frac(:,:,j) = 1.0
656	ENDIF
657	ENDDO
658	ENDIF
659	!end gmjc
660	!! 2.2.6 for security: can't add more lignin than structural litter (above/below)
661	DO l = 1, nlevs !Loop over litter levels (above and below ground)
662	DO m = 1,nvm !Lopp over PFTs
663
664	lignin_struc_inc(:,m,l) = &
665	MIN( lignin_struc_inc(:,m,l), litter_inc(:,istructural,m,l,icarbon) )
666
667	ENDDO
668	ENDDO
669
670	!! 2.2.7 new lignin content: add old lignin and lignin increase, divide by
671	!! total structural litter (above/below)
672	DO l = 1, nlevs !Loop over litter levels (above and below ground)
673	DO m = 1,nvm !Loop over PFTs
674	WHERE( litter(:,istructural,m,l,icarbon) .GT. min_stomate )
675
676	!MM : Soenke modif
677	! Best vectorization ?
678	!!$ lignin_struc(:,:,:) = &
679	!!$ ( lignin_struc(:,:,:)*old_struc(:,:,:) + lignin_struc_inc(:,:,:) ) / &
680	!!$ litter(:,istructural,:,:,icarbon)
681
682	lignin_struc(:,m,l) = lignin_struc(:,m,l) * old_struc(:,m,l)
683	lignin_struc(:,m,l) = lignin_struc(:,m,l) + lignin_struc_inc(:,m,l)
684	lignin_struc(:,m,l) = lignin_struc(:,m,l) / litter(:,istructural,m,l,icarbon)
685	ELSEWHERE
686	lignin_struc(:,m,l) = zero
687	ENDWHERE
688	ENDDO
689	ENDDO
690
691
692	!! 2.3 new litter fraction per PFT (for structural and metabolic litter, above
693	!! the ground).
694	DO j = 1,nvm !Loop over PFTs
695
696	WHERE ( litter(:,:,j,iabove,icarbon) .GT. min_stomate )
697
698	litterpart(:,j,:) = &
699	( litter_pft(:,j,:,icarbon) + litter_inc_PFT(:,j,:,iabove,icarbon) ) / litter(:,:,j,iabove,icarbon)
700
701	ELSEWHERE
702
703	litterpart(:,j,:) = zero
704
705	ENDWHERE
706
707	ENDDO
708
709
710	!! 3. Temperature control on decay: Factor between 0 and 1
711
712	!! 3.1 above: surface temperature
713	control_temp(:,iabove) = control_temp_func (npts,tsurf,frozen_respiration_func)
714
715	!!!qcj++ peatland: in control_temp_func_peat, changed Q10 value for peat vegetation
716	IF (perma_peat) THEN
717	control_temp_peat(:,iabove)=control_temp_func_peat (npts,tsurf,frozen_respiration_func)
718	ENDIF
719	!! 3.2 below: convolution of temperature and decomposer profiles
720	!! (exponential decomposer profile supposed)
721
722	!! 3.2.1 rpc is an integration constant such that the integral of the root profile is 1.
723	rpc(:) = un / ( un - EXP( -z_soil(nslm) / z_decomp ) )
724
725	!! 3.2.2 integrate over the nbdl levels
726	tsoil_decomp(:) = zero
727
728	DO l = 1, nslm
729
730	tsoil_decomp(:) = &
731	tsoil_decomp(:) + tsoil(:,l) * rpc(:) * &
732	( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
733
734	ENDDO
735
736	control_temp(:,ibelow) = control_temp_func (npts,tsoil_decomp,frozen_respiration_func)
737
738	!!!qcj++ peatland: in control_temp_func_peat, changed Q10 value for peat vegetation
739	IF (perma_peat) THEN
740	control_temp_peat(:,ibelow)=control_temp_func_peat (npts,tsoil_decomp,frozen_respiration_func)
741	ENDIF
742
743
744	!!!! crops
745
746	! for STICS
747	! 3.3 calculation of soil temperature with the resolution of 1 cm
748	! here, we calculated the soil temperature with 1 cm for three layers around sowing depth, see STICS code
749
750	! initialization of the tsoil_cm
751
752	tsoil_cm(:, :, :) = zero
753	ldep(:, :) = zero ! different for different PFTs
754	! the depths of three layers around sowing depth, in cm and they are transformed into integer
755
756	DO ipft = 2, nvm
757	IF (ok_LAIdev(ipft)) THEN
758	ldep(ipft, :) = (/SP_profsem(ipft) - 1.0, SP_profsem(ipft), SP_profsem(ipft) + 1.0/) !three depths, SP_profsem is in cma
759	!!! 3cm by default
760	ENDIF
761	ENDDO
762
763	IF (printlev>=4) THEN
764	WRITE(,) 'three layer depth in litter is', ldep(13, :)
765	WRITE(,) 'z_soil in litter is', z_soil
766	ENDIF
767	! loop for each layer of the three layers and calculate the soil temperature at the specific depth
768
769
770	DO ipft = 2, nvm
771	IF (ok_LAIdev(ipft)) THEN
772
773	DO il = 1, SIZE(ldep, 2)! loop for the three layers, three depth
774
775	! initialize the surface soil depth, that means 0 cm
776	uldep = z_soil(0)*100.0 ! the upper limit of the first layer, and also the initial depth
777	bldep = z_soil(1)/2.0*100.0 ! the lower limit of the first half layer, and also the initial depth
778	! initialize the surface soil temperature
779	DO nhl = 1, 2*nslm
780	ultsoil(:, nhl) = tsurf(:) ! the surface soil temperature and also the initial soil temperature
781	ENDDO
782
783	DO nhl = 1, 2*nslm
784	IF ((ldep(ipft, il) .GE. uldep) .AND. (ldep(ipft, il) .LT. bldep)) THEN !
785	IF (nhl < 2*nslm) THEN
786	tsoil_cm(:,ipft, il) = &
787	& ultsoil(:,nhl) + &
788	& ((tsoil(:, nhl - FLOOR((nhl-1.0)/2.0)) + tsoil(:, NINT(nhl/2.0)))/2.0 - ultsoil(:, nhl))&
789	& /(bldep - uldep)&
790	& *(ldep(ipft, il) - uldep)
791	! linear interpolation
792	ELSE
793	tsoil_cm(:, ipft, il) = tsoil(:, nslm)
794	ENDIF
795	ENDIF
796
797	uldep = bldep
798	bldep = bldep + (z_soil(CEILING(nhl/2.0))100.0 - z_soil(CEILING(nhl/2.0)-1)100.0)/2.0
799	IF (nhl < 2*nslm) THEN
800	ultsoil(:, nhl) =(tsoil(:, nhl - FLOOR((nhl-1.0)/2.0)) + tsoil(:, NINT(nhl/2.0)))/2.0 !
801	ELSE
802	ultsoil(:, nhl) = tsoil(:, nslm)
803	ENDIF
804	ENDDO ! nhl
805	ENDDO ! il
806	ENDIF
807	ENDDO ! ipft
808
809
810	!!!! end crops, xuhui
811
812	!! 4. Moisture control. Factor between 0 and 1
813
814	!! 4.1 above the ground: litter humidity
815	control_moist(:,iabove) = control_moist_func (npts, litterhum)
816	IF (perma_peat) THEN
817
818	!!!qcj++ peatland : use water content instead of relative hum
819	control_moist_peat(:,iabove)=control_moist_func_peat(npts, mc_peat_above)
820	ENDIF
821	!
822	!! 4.2 below: convolution of humidity and decomposer profiles
823	! (exponential decomposer profile supposed)
824
825	!! 4.2.1 rpc is an integration constant such that the integral of the root profile is 1.
826	rpc(:) = un / ( un - EXP( -z_soil(nslm) / z_decomp ) )
827
828	!! 4.2.2 integrate over the nbdl levels
829	soilhum_decomp(:) = zero
830
831	DO l = 1, nslm !Loop over soil levels
832
833	soilhum_decomp(:) = &
834	soilhum_decomp(:) + soilhum(:,l) * rpc(:) * &
835	( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
836
837	ENDDO
838
839	control_moist(:,ibelow) = control_moist_func (npts, soilhum_decomp)
840
841	!!!qcj++ peatland : use water content instead of relative hum
842	IF (perma_peat) THEN
843	mc_peat_below(:) =zero
844	DO l = 1, nslm !Loop over soil levels
845	mc_peat_below(:) = &
846	mc_peat_below(:)+shumdiag_peat(:,l) * rpc(:) * &
847	( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
848	ENDDO
849	control_moist_peat(:,ibelow)=control_moist_func_peat(npts, mc_peat_below)
850	ENDIF
851	!!!!! crops
852
853	! 4.3 calculation of soil relative humidity of three depths
854	! similar to the calculation of soil temperature, see above
855
856	! initialize the soil humidity for the three layers
857
858	soilhum_cm(:, :, :) = zero
859	ldep(:, :) = zero
860
861
862	! the depths for the three layers around the sowing depth
863
864	DO ipft = 1, nvm
865	IF (ok_LAIdev(ipft)) THEN
866	ldep(ipft, :) = (/SP_profsem(ipft) - 1.0, SP_profsem(ipft), SP_profsem(ipft) + 1.0/) !three depths, SP_profsem is in cm
867	ENDIF
868	ENDDO
869
870	! IF (printlev>=4) WRITE(,) 'in stomate_litter, the ldep of wheat is', ldep(12, :)
871
872
873	! loop calculation
874
875	DO ipft = 1, nvm
876	IF (ok_LAIdev(ipft)) THEN
877	DO il = 1, SIZE(ldep, 2)!the ith layer
878
879	! initialize the surface depths, this variable can be changed for reassign a value for the upper limit depth for the ith half layer
880	!uldep(:) = z_soil(1)*100.0 ! the upper limit of the first layer, and also the initial depth
881
882	uldep = z_soil(0)*100.0 ! the upper limit of the first layer, and also the initial depth
883	bldep = z_soil(1)/2.0*100.0
884
885	! initialize the surface soil moisture, is it correct?
886	DO nhl = 1, 2*nslm
887	ulhum(:, nhl) = soilhum(:, 1) ! the surface soil moisture and also the initial soil moisture
888	ENDDO
889
890	DO nhl = 1, 2*nslm
891	IF ((ldep(ipft, il) .GE. uldep) .AND. (ldep(ipft,il) .LT. bldep)) THEN
892	IF (nhl < 2*nslm) THEN
893	soilhum_cm(:, ipft, il) = &
894	& ulhum(:, nhl) + &
895	& ((soilhum(:, nhl - FLOOR((nhl-1)/2.0)) + soilhum(:, NINT(nhl/2.0)))/2.0 - ulhum(:, nhl))&
896	& /(bldep - uldep)&
897	& *(ldep(ipft, il) - uldep)
898	! linear interpolation
899	ELSE
900	soilhum_cm(:, ipft, il) = soilhum(:, nslm)
901	ENDIF
902	ENDIF
903
904	uldep = bldep
905	bldep = bldep + (z_soil(CEILING(nhl/2.0))100.0 - z_soil(CEILING(nhl/2.0)-1)100.0)/2.0
906	IF (nhl < 2*nslm) THEN
907	ulhum(:, nhl) =(soilhum(:, nhl - FLOOR((nhl-1.0)/2.0)) + soilhum(:, NINT(nhl/2.0)))/2.0 !
908	ELSE
909	ulhum(:, nhl) = soilhum(:, nslm)
910	ENDIF
911	ENDDO ! nhl
912	ENDDO ! il
913	ENDIF
914	ENDDO !ipft
915	!!!!! end crops, xuhui
916
917	!! 5. fluxes from litter to carbon pools and respiration
918
919	DO l = 1, nlevs !Loop over litter levels (above and below ground)
920	DO m = 1,nvm !Loop over PFTs
921
922	!! 5.1 structural litter: goes into active and slow carbon pools + respiration
923
924	!! 5.1.1 total quantity of structural litter which is decomposed
925	!!!qcj++ peatland
926	IF (perma_peat .AND. is_peat(m)) THEN
927	fd(:) = dt/litter_tau(istructural) * &
928	control_temp_peat(:,l) * control_moist_peat(:,l) * exp(-litter_struct_coef * lignin_struc(:,m,l) )
929	ELSE
930	fd(:) = dt/litter_tau(istructural) * &
931	control_temp(:,l) * control_moist(:,l) * exp( -litter_struct_coef * lignin_struc(:,m,l) )
932	ENDIF
933
934	DO k = 1,nelements
935
936	qd(:,k) = litter(:,istructural,m,l,k) * fd(:)
937
938	END DO
939
940	litter(:,istructural,m,l,:) = litter(:,istructural,m,l,:) - qd(:,:)
941	!spitfire
942	!update the fuel poop size accordingly after respiration
943	qd_fuel(:,:)=qd(:,:)
944
945	IF(l==iabove) THEN
946	DO k=1, npts
947	IF (fuel_nlitt_total_pft(k,m,istructural,icarbon).GT.min_stomate) THEN
948
949	fuel_1hr(k,m,istructural,icarbon) = fuel_1hr(k,m,istructural,icarbon) - (qd_fuel(k,icarbon) * &
950	fuel_1hr(k,m,istructural,icarbon)/fuel_nlitt_total_pft(k,m,istructural,icarbon))
951
952	fuel_10hr(k,m,istructural,icarbon) = fuel_10hr(k,m,istructural,icarbon) - (qd_fuel(k,icarbon) * &
953	fuel_10hr(k,m,istructural,icarbon)/fuel_nlitt_total_pft(k,m,istructural,icarbon))
954
955	fuel_100hr(k,m,istructural,icarbon) = fuel_100hr(k,m,istructural,icarbon) - (qd_fuel(k,icarbon) * &
956	fuel_100hr(k,m,istructural,icarbon)/fuel_nlitt_total_pft(k,m,istructural,icarbon))
957
958	fuel_1000hr(k,m,istructural,icarbon) = fuel_1000hr(k,m,istructural,icarbon) - (qd_fuel(k,icarbon) * &
959	fuel_1000hr(k,m,istructural,icarbon)/fuel_nlitt_total_pft(k,m,istructural,icarbon))
960
961	fuel_nlitt_total_pft(k,m,istructural,icarbon) = fuel_1hr(k,m,istructural,icarbon) + &
962	fuel_10hr(k,m,istructural,icarbon) + fuel_100hr(k,m,istructural,icarbon) + &
963	fuel_1000hr(k,m,istructural,icarbon)
964	ENDIF
965
966	diff_frac(k) = litter(k,istructural,m,l,icarbon) - fuel_nlitt_total_pft(k,m,istructural,icarbon)
967
968	IF (ABS(diff_frac(k)).GT.min_stomate .AND. fuel_nlitt_total_pft(k,m,istructural,icarbon).GT.min_stomate) THEN
969	! simply divided even fraction to each fuel class, that is 1/4 of buffer
970	fuel_1hr(k,m,istructural,icarbon) = fuel_1hr(k,m,istructural,icarbon) + fuel_1hr(k,m,istructural,icarbon) &
971	* diff_frac(k)/fuel_nlitt_total_pft(k,m,istructural,icarbon)
972	fuel_10hr(k,m,istructural,icarbon) = fuel_10hr(k,m,istructural,icarbon) + fuel_10hr(k,m,istructural,icarbon) &
973	* diff_frac(k)/fuel_nlitt_total_pft(k,m,istructural,icarbon)
974	fuel_100hr(k,m,istructural,icarbon) = fuel_100hr(k,m,istructural,icarbon) + fuel_100hr(k,m,istructural,icarbon) &
975	* diff_frac(k)/fuel_nlitt_total_pft(k,m,istructural,icarbon)
976	fuel_1000hr(k,m,istructural,icarbon) = fuel_1000hr(k,m,istructural,icarbon) + fuel_1000hr(k,m,istructural,icarbon) &
977	* diff_frac(k)/fuel_nlitt_total_pft(k,m,istructural,icarbon)
978	fuel_nlitt_total_pft(k,m,istructural,icarbon) = fuel_1hr(k,m,istructural,icarbon) + &
979	fuel_10hr(k,m,istructural,icarbon) + fuel_100hr(k,m,istructural,icarbon) + &
980	fuel_1000hr(k,m,istructural,icarbon)
981	ENDIF
982	ENDDO
983	#ifdef STRICT_CHECK
984	! Check consistency in fuel_Xhr variables. Negative values are not allowed
985	IF (ANY(fuel_1hr < 0)) CALL ipslerr_p(3, 'stomate_litter', 'fuel_1hr has negative values', '', '')
986	IF (ANY(fuel_10hr < 0)) CALL ipslerr_p(3, 'stomate_litter', 'fuel_10hr has negative values', '', '')
987	IF (ANY(fuel_100hr < 0)) CALL ipslerr_p(3, 'stomate_litter', 'fuel_100hr has negative values', '', '')
988	IF (ANY(fuel_1000hr < 0)) CALL ipslerr_p(3, 'stomate_litter', 'fuel_1000hr has negative values', '', '')
989	#endif
990	ENDIF
991
992	!endspit
993
994	!! 5.1.2 decompose same fraction of structural part of dead leaves. Not exact
995	!! as lignine content is not the same as that of the total structural litter.
996	! to avoid a multiple (for ibelow and iabove) modification of dead_leaves,
997	! we do this test to do this calcul only ones in 1,nlev loop
998	if (l == iabove) dead_leaves(:,m,istructural) = dead_leaves(:,m,istructural) * ( un - fd(:) )
999
1000	!! 5.1.3 non-lignin fraction of structural litter goes into
1001	!! active carbon pool + respiration
1002	soilcarbon_input(:,iactive,m) = soilcarbon_input(:,iactive,m) + &
1003	frac_soil(istructural,iactive,l) * qd(:,icarbon) * ( 1. - lignin_struc(:,m,l) ) / dt
1004	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
1005	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
1006	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dt_sechiba/one_day
1007	! time step and avoid some constructions that could create bug during future developments.
1008	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
1009	( 1. - frac_soil(istructural,iactive,l) ) * qd(:,icarbon) * &
1010	( 1. - lignin_struc(:,m,l) ) / dt
1011
1012	!! 5.1.4 lignin fraction of structural litter goes into
1013	!! slow carbon pool + respiration
1014	soilcarbon_input(:,islow,m) = soilcarbon_input(:,islow,m) + &
1015	frac_soil(istructural,islow,l) * qd(:,icarbon) * lignin_struc(:,m,l) / dt
1016
1017	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
1018	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
1019	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dt_sechiba/one_day
1020	! time step and avoid some constructions that could create bug during future developments.
1021	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
1022	( 1. - frac_soil(istructural,islow,l) ) * qd(:,icarbon) * lignin_struc(:,m,l) / dt
1023
1024
1025	!! 5.2 metabolic litter goes into active carbon pool + respiration
1026
1027	!! 5.2.1 total quantity of metabolic litter that is decomposed
1028	!!!qcj++ peatland
1029	IF (perma_peat .AND. is_peat(m)) THEN
1030	fd(:) = dt/litter_tau(imetabolic) * control_temp_peat(:,l) * control_moist_peat(:,l)
1031	ELSE
1032	fd(:) = dt/litter_tau(imetabolic) * control_temp(:,l) * control_moist(:,l)
1033	ENDIF
1034
1035	DO k = 1,nelements
1036
1037	qd(:,k) = litter(:,imetabolic,m,l,k) * fd(:)
1038
1039	END DO
1040
1041	litter(:,imetabolic,m,l,:) = litter(:,imetabolic,m,l,:) - qd(:,:)
1042
1043	!spitfire
1044	!update aboveground fuel load
1045	qd_fuel(:,:) = qd(:,:)
1046
1047	IF(l==iabove) THEN
1048	DO k=1, npts
1049	IF (fuel_nlitt_total_pft(k,m,imetabolic,icarbon).gt.min_stomate) THEN
1050	fuel_1hr(k,m,imetabolic,icarbon) = fuel_1hr(k,m,imetabolic,icarbon) - (qd_fuel(k,icarbon) * &
1051	fuel_1hr(k,m,imetabolic,icarbon)/fuel_nlitt_total_pft(k,m,imetabolic,icarbon))
1052
1053	fuel_10hr(k,m,imetabolic,icarbon) = fuel_10hr(k,m,imetabolic,icarbon) - (qd_fuel(k,icarbon) * &
1054	fuel_10hr(k,m,imetabolic,icarbon)/fuel_nlitt_total_pft(k,m,imetabolic,icarbon))
1055
1056	fuel_100hr(k,m,imetabolic,icarbon) = fuel_100hr(k,m,imetabolic,icarbon) - (qd_fuel(k,icarbon) * &
1057	fuel_100hr(k,m,imetabolic,icarbon)/fuel_nlitt_total_pft(k,m,imetabolic,icarbon))
1058
1059	fuel_1000hr(k,m,imetabolic,icarbon) = fuel_1000hr(k,m,imetabolic,icarbon) - (qd_fuel(k,icarbon) * &
1060	fuel_1000hr(k,m,imetabolic,icarbon)/fuel_nlitt_total_pft(k,m,imetabolic,icarbon))
1061
1062	fuel_nlitt_total_pft(k,m,imetabolic,icarbon) = fuel_1hr(k,m,imetabolic,icarbon) + &
1063	fuel_10hr(k,m,imetabolic,icarbon) + fuel_100hr(k,m,imetabolic,icarbon) + &
1064	fuel_1000hr(k,m,imetabolic,icarbon)
1065	ENDIF
1066
1067	diff_frac(k) = litter(k,imetabolic,m,l,icarbon) - fuel_nlitt_total_pft(k,m,imetabolic,icarbon)
1068
1069	IF(ABS(diff_frac(k)).GT.min_stomate .AND. fuel_nlitt_total_pft(k,m,imetabolic,icarbon).gt.min_stomate) THEN
1070	! simply divided even fraction to each fuel class, that is 1/4 of buffer
1071	fuel_1hr(k,m,imetabolic,icarbon) = fuel_1hr(k,m,imetabolic,icarbon) + fuel_1hr(k,m,imetabolic,icarbon)&
1072	* diff_frac(k)/fuel_nlitt_total_pft(k,m,imetabolic,icarbon)
1073	fuel_10hr(k,m,imetabolic,icarbon) = fuel_10hr(k,m,imetabolic,icarbon) + fuel_10hr(k,m,imetabolic,icarbon)&
1074	* diff_frac(k)/fuel_nlitt_total_pft(k,m,imetabolic,icarbon)
1075	fuel_100hr(k,m,imetabolic,icarbon) = fuel_100hr(k,m,imetabolic,icarbon) + fuel_100hr(k,m,imetabolic,icarbon)&
1076	* diff_frac(k)/fuel_nlitt_total_pft(k,m,imetabolic,icarbon)
1077	fuel_1000hr(k,m,imetabolic,icarbon) = fuel_1000hr(k,m,imetabolic,icarbon) + fuel_1000hr(k,m,imetabolic,icarbon)&
1078	* diff_frac(k)/fuel_nlitt_total_pft(k,m,imetabolic,icarbon)
1079	fuel_nlitt_total_pft(k,m,imetabolic,icarbon) = fuel_1hr(k,m,imetabolic,icarbon) + &
1080	fuel_10hr(k,m,imetabolic,icarbon) + fuel_100hr(k,m,imetabolic,icarbon) + &
1081	fuel_1000hr(k,m,imetabolic,icarbon)
1082	ENDIF
1083	ENDDO
1084	#ifdef STRICT_CHECK
1085	! Check consistency in fuel_Xhr variables. Negative values are not allowed
1086	IF (ANY(fuel_1hr < 0)) CALL ipslerr_p(3, 'stomate_litter', 'fuel_1hr has negative values', '', '')
1087	IF (ANY(fuel_10hr < 0)) CALL ipslerr_p(3, 'stomate_litter', 'fuel_10hr has negative values', '', '')
1088	IF (ANY(fuel_100hr < 0)) CALL ipslerr_p(3, 'stomate_litter', 'fuel_100hr has negative values', '', '')
1089	IF (ANY(fuel_1000hr < 0)) CALL ipslerr_p(3, 'stomate_litter', 'fuel_1000hr has negative values', '', '')
1090	#endif
1091	ENDIF
1092	!endspit
1093
1094	!! 5.2.2 decompose same fraction of metabolic part of dead leaves.
1095	! to avoid a multiple (for ibelow and iabove) modification of dead_leaves,
1096	! we do this test to do this calcul only ones in 1,nlev loop
1097	if (l == iabove) dead_leaves(:,m,imetabolic) = dead_leaves(:,m,imetabolic) * ( 1. - fd(:) )
1098
1099
1100	!! 5.2.3 put decomposed litter into carbon pool + respiration
1101	soilcarbon_input(:,iactive,m) = soilcarbon_input(:,iactive,m) + &
1102	frac_soil(imetabolic,iactive,l) * qd(:,icarbon) / dt
1103
1104	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
1105	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
1106	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dt_sechiba/one_day
1107	! time step and avoid some constructions that could create bug during future developments.
1108	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
1109	( 1. - frac_soil(imetabolic,iactive,l) ) * qd(:,icarbon) / dt
1110
1111	ENDDO
1112	ENDDO
1113
1114
1115	!! 6. calculate fraction of total soil covered by dead leaves
1116
1117	CALL deadleaf (npts, veget_max, dead_leaves, deadleaf_cover, &!)
1118	!gmjc
1119	sla_calc)
1120	!end gmjc
1121
1122	!! 7. (Quasi-)Analytical Spin-up : Start filling MatrixA
1123
1124	IF (spinup_analytic) THEN
1125
1126	MatrixA(:,:,:,:) = zero
1127	VectorB(:,:,:) = zero
1128
1129
1130	DO m = 1,nvm
1131
1132	!- MatrixA : carbon fluxes leaving the litter
1133
1134	MatrixA(:,m,istructural_above,istructural_above)= - dt/litter_tau(istructural) * &
1135	control_temp(:,iabove) * control_moist(:,iabove) * exp( -litter_struct_coef * lignin_struc(:,m,iabove) )
1136
1137	MatrixA(:,m,istructural_below,istructural_below) = - dt/litter_tau(istructural) * &
1138	control_temp(:,ibelow) * control_moist(:,ibelow) * exp( -litter_struct_coef * lignin_struc(:,m,ibelow) )
1139
1140	MatrixA(:,m,imetabolic_above,imetabolic_above) = - dt/litter_tau(imetabolic) * &
1141	control_temp(:,iabove) * control_moist(:,iabove)
1142
1143	MatrixA(:,m,imetabolic_below,imetabolic_below) = - dt/litter_tau(imetabolic) * &
1144	control_temp(:,ibelow) * control_moist(:,ibelow)
1145
1146
1147	!- MatrixA : carbon fluxes between the litter and the pools (the rest of the matrix is filled in stomate_soilcarbon.f90)
1148
1149	MatrixA(:,m,iactive_pool,istructural_above) = frac_soil(istructural,iactive,iabove) * &
1150	dt/litter_tau(istructural) * &
1151	control_temp(:,iabove) * control_moist(:,iabove) * &
1152	exp( -litter_struct_coef * lignin_struc(:,m,iabove) ) * &
1153	( 1. - lignin_struc(:,m,iabove) )
1154
1155
1156	MatrixA(:,m,iactive_pool,istructural_below) = frac_soil(istructural,iactive,ibelow) * &
1157	dt/litter_tau(istructural) * &
1158	control_temp(:,ibelow) * control_moist(:,ibelow) * &
1159	exp( -litter_struct_coef * lignin_struc(:,m,ibelow) ) * &
1160	( 1. - lignin_struc(:,m,ibelow) )
1161
1162
1163	MatrixA(:,m,iactive_pool,imetabolic_above) = frac_soil(imetabolic,iactive,iabove) * &
1164	dt/litter_tau(imetabolic) * control_temp(:,iabove) * control_moist(:,iabove)
1165
1166
1167	MatrixA(:,m,iactive_pool,imetabolic_below) = frac_soil(imetabolic,iactive,ibelow) * &
1168	dt/litter_tau(imetabolic) * control_temp(:,ibelow) * control_moist(:,ibelow)
1169
1170	MatrixA(:,m,islow_pool,istructural_above) = frac_soil(istructural,islow,iabove) * &
1171	dt/litter_tau(istructural) * &
1172	control_temp(:,iabove) * control_moist(:,iabove) * &
1173	exp( -litter_struct_coef * lignin_struc(:,m,iabove) )* &
1174	lignin_struc(:,m,iabove)
1175
1176
1177	MatrixA(:,m,islow_pool,istructural_below) = frac_soil(istructural,islow,ibelow) * &
1178	dt/litter_tau(istructural) * &
1179	control_temp(:,ibelow) * control_moist(:,ibelow) * &
1180	exp( -litter_struct_coef * lignin_struc(:,m,ibelow) )* &
1181	lignin_struc(:,m,ibelow)
1182
1183
1184	!- VectorB : carbon input -
1185
1186	VectorB(:,m,istructural_above) = litter_inc_PFT(:,m,istructural,iabove,icarbon)
1187	VectorB(:,m,istructural_below) = litter_inc_PFT(:,m,istructural,ibelow,icarbon)
1188	VectorB(:,m,imetabolic_above) = litter_inc_PFT(:,m,imetabolic,iabove,icarbon)
1189	VectorB(:,m,imetabolic_below) = litter_inc_PFT(:,m,imetabolic,ibelow,icarbon)
1190
1191	IF (printlev>=4) WRITE(numout,*) 'We filled MatrixA and VectorB'
1192
1193	ENDDO ! Loop over # PFTs
1194
1195	ENDIF ! spinup analytic
1196	IF (printlev>=4) WRITE(numout,*) 'Leaving littercalc'
1197
1198	END SUBROUTINE littercalc
1199
1200
1201	!! ==============================================================================================================================\n
1202	!! SUBROUTINE : deadleaf
1203	!!
1204	!>\BRIEF This routine calculates the deadleafcover.
1205	!!
1206	!! DESCRIPTION : It first calculates the lai corresponding to the dead leaves (LAI) using
1207	!! the dead leaves carbon content (DL) the specific leaf area (sla) and the
1208	!! maximal coverage fraction of a PFT (vegetmax) using the following equations:
1209	!! \latexonly
1210	!! \input{deadleaf1.tex}
1211	!! \endlatexonly
1212	!! \n
1213	!! Then, the dead leaf cover (DLC) is calculated as following:\n
1214	!! \latexonly
1215	!! \input{deadleaf2.tex}
1216	!! \endlatexonly
1217	!! \n
1218	!!
1219	!! RECENT CHANGE(S) : None
1220	!!
1221	!! MAIN OUTPUT VARIABLE: ::deadleaf_cover
1222	!!
1223	!! REFERENCE(S) : None
1224	!!
1225	!! FLOWCHART : None
1226	!! \n
1227	!_ ================================================================================================================================
1228
1229	SUBROUTINE deadleaf (npts, veget_max, dead_leaves, deadleaf_cover, &
1230	!gmjc
1231	sla_calc)
1232	!end gmjc
1233	!! 0. Variable and parameter declaration
1234
1235	!! 0.1 Input variables
1236
1237	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixels (unitless)
1238	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: dead_leaves !! Dead leaves per ground unit area, per PFT,
1239	!! metabolic and structural
1240	!! @tex $(gC m^{-2})$ @endtex
1241	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT defined in
1242	!! the input vegetation map
1243	!! @tex $(m^2 m^{-2})$ @endtex
1244
1245	!! 0.2 Output variables
1246
1247	REAL(r_std), DIMENSION(npts), INTENT(out) :: deadleaf_cover !! Fraction of soil covered by dead leaves over all PFTs
1248	!! (0-1, unitless)
1249	!gmjc
1250	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: sla_calc
1251	!end gmjc
1252	!! 0.3 Modified variables
1253
1254	!! 0.4 Local variables
1255
1256	REAL(r_std), DIMENSION(npts) :: dead_lai !! LAI of dead leaves @tex $(m^2 m^{-2})$ @endtex
1257	INTEGER(i_std) :: j !! Index (unitless)
1258	!_ ================================================================================================================================
1259
1260	!! 1. LAI of dead leaves
1261
1262	dead_lai(:) = zero
1263
1264	DO j = 1,nvm !Loop over PFTs
1265	!JCMODIF
1266	! dead_lai(:) = dead_lai(:) + ( dead_leaves(:,j,imetabolic) + dead_leaves(:,j,istructural) ) * sla(j) &
1267	dead_lai(:) = dead_lai(:) + ( dead_leaves(:,j,imetabolic) + &
1268	dead_leaves(:,j,istructural) ) * sla_calc(:,j) &
1269	!ENDJCMODIF
1270	* veget_max(:,j)
1271	ENDDO
1272
1273	!! 2. fraction of soil covered by dead leaves
1274
1275	deadleaf_cover(:) = un - exp( - 0.5 * dead_lai(:) )
1276
1277	IF (printlev>=4) WRITE(numout,*) 'Leaving deadleaf'
1278
1279	END SUBROUTINE deadleaf
1280
1281
1282	!! ================================================================================================================================
1283	!! FUNCTION : control_moist_func
1284	!!
1285	!>\BRIEF Calculate moisture control for litter and soil C decomposition
1286	!!
1287	!! DESCRIPTION : Calculate moisture control factor applied
1288	!! to litter decomposition and to soil carbon decomposition in
1289	!! stomate_soilcarbon.f90 using the following equation: \n
1290	!! \latexonly
1291	!! \input{control_moist_func1.tex}
1292	!! \endlatexonly
1293	!! \n
1294	!! with M the moisture control factor and soilmoisutre, the soil moisture
1295	!! calculated in sechiba.
1296	!! Then, the function is ranged between Moistcont_min and 1:\n
1297	!! \latexonly
1298	!! \input{control_moist_func2.tex}
1299	!! \endlatexonly
1300	!! \n
1301	!! RECENT CHANGE(S) : None
1302	!!
1303	!! RETURN VALUE : ::moistfunc_result
1304	!!
1305	!! REFERENCE(S) : None
1306	!!
1307	!! FLOWCHART : None
1308	!! \n
1309	!_ ================================================================================================================================
1310
1311	FUNCTION control_moist_func (npts, moist_in) RESULT (moistfunc_result)
1312
1313	!! 0. Variable and parameter declaration
1314
1315	!! 0.1 Input variables
1316
1317	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixel (unitless)
1318	REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_in !! relative humidity (unitless)
1319
1320	!! 0.2 Output variables
1321
1322	REAL(r_std), DIMENSION(npts) :: moistfunc_result !! Moisture control factor (0.25-1, unitless)
1323
1324	!! 0.3 Modified variables
1325
1326	!! 0.4 Local variables
1327
1328	!_ ================================================================================================================================
1329
1330	moistfunc_result(:) = -moist_coeff(1) * moist_in(:) * moist_in(:) + moist_coeff(2)* moist_in(:) - moist_coeff(3)
1331	moistfunc_result(:) = MAX( moistcont_min, MIN( un, moistfunc_result(:) ) )
1332
1333	END FUNCTION control_moist_func
1334
1335
1336	!! ================================================================================================================================
1337	!! FUNCTION : control_temp_func
1338	!!
1339	!>\BRIEF Calculate temperature control for litter and soild C decomposition
1340	!!
1341	!! DESCRIPTION : Calculate temperature control factor applied
1342	!! to litter decomposition and to soil carbon decomposition in
1343	!! stomate_soilcarbon.f90 using the following equation: \n
1344	!! \latexonly
1345	!! \input{control_temp_func1.tex}
1346	!! \endlatexonly
1347	!! \n
1348	!! with T the temperature control factor, temp the temperature in Kelvin of
1349	!! the air (for aboveground litter) or of the soil (for belowground litter
1350	!! and soil)
1351	!! Then, the function is limited in its maximal range to 1:\n
1352	!! \latexonly
1353	!! \input{control_temp_func2.tex}
1354	!! \endlatexonly
1355	!! \n
1356	!! RECENT CHANGE(S) : None
1357	!!
1358	!! RETURN VALUE: ::tempfunc_result
1359	!!
1360	!! REFERENCE(S) : None
1361	!!
1362	!! FLOWCHART : None
1363	!! \n
1364	!_ ================================================================================================================================
1365
1366	FUNCTION control_temp_func (npts, temp_in, frozen_respiration_func) RESULT (tempfunc_result)
1367
1368	!! 0. Variable and parameter declaration
1369
1370	!! 0.1 Input variables
1371	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of land pixels (unitless)
1372	REAL(r_std), DIMENSION(npts), INTENT(in) :: temp_in !! Temperature (K)
1373	INTEGER(i_std), INTENT(in) :: frozen_respiration_func
1374	!! 0.2 Output variables
1375	REAL(r_std), DIMENSION(npts) :: tempfunc_result !! Temperature control factor (0-1, unitless)
1376
1377	!! 0.3 Modified variables
1378
1379	!! 0.4 Local variables
1380
1381	!_ ================================================================================================================================
1382	SELECT CASE(frozen_respiration_func)
1383
1384	CASE(0) !!! this is the standard ORCHIDEE state
1385
1386	tempfunc_result(:) = exp( soil_Q10 * ( temp_in(:) - (ZeroCelsius+tsoil_ref)) / Q10 )
1387	tempfunc_result(:) = MIN( un, tempfunc_result(:) )
1388
1389	CASE(1) !!! cutoff respiration when T < -1C
1390	WHERE (temp_in(:) .GT. ZeroCelsius ) !!! normal as above
1391	tempfunc_result(:) = EXP( 0.69 * ( temp_in(:) - (ZeroCelsius+30.) ) / 10. )
1392	ELSEWHERE (temp_in(:) .GT. ZeroCelsius - 1. ) !!! linear dropoff to zero
1393	tempfunc_result(:) = (temp_in(:) - (ZeroCelsius - 1.)) * EXP( 0.69 * ( ZeroCelsius - (ZeroCelsius+30.) ) / 10. )
1394	ELSEWHERE !!! zero
1395	tempfunc_result(:) = zero
1396	endwhere
1397
1398	tempfunc_result(:) = MAX(MIN( 1._r_std, tempfunc_result(:) ), zero)
1399
1400	CASE(2) !!! cutoff respiration when T < -3C
1401	WHERE (temp_in(:) .GT. ZeroCelsius ) !!! normal as above
1402	tempfunc_result(:) = EXP( 0.69 * ( temp_in(:) - (ZeroCelsius+30.) ) / 10. )
1403	ELSEWHERE (temp_in(:) .GT. ZeroCelsius - 3. ) !!! linear dropoff to zero
1404	tempfunc_result(:) = ((temp_in(:) - (ZeroCelsius - 3.))/3.) * EXP( 0.69 * ( ZeroCelsius - (ZeroCelsius+30.) ) / 10. )
1405	ELSEWHERE !!! zero
1406	tempfunc_result(:) = zero
1407	endwhere
1408
1409	CASE(3) !!! q10 = 100 when below zero
1410	WHERE (temp_in(:) .GT. ZeroCelsius ) !!! normal as above
1411	tempfunc_result(:) = EXP( 0.69 * ( temp_in(:) - (ZeroCelsius+30.) ) / 10. )
1412	ELSEWHERE
1413	tempfunc_result(:) = EXP( 4.605 * ( temp_in(:) - (ZeroCelsius) ) / 10.) * EXP( 0.69 * ( -30. ) / 10. )
1414	endwhere
1415
1416	CASE(4) !!! q10 = 1000 when below zero
1417	WHERE (temp_in(:) .GT. ZeroCelsius ) !!! normal as above
1418	tempfunc_result(:) = EXP( 0.69 * ( temp_in(:) - (ZeroCelsius+30.) ) / 10. )
1419	ELSEWHERE
1420	tempfunc_result(:) = EXP( 6.908 * ( temp_in(:) - (ZeroCelsius) ) / 10.) * EXP( 0.69 * ( -30. ) / 10. )
1421	endwhere
1422
1423	END SELECT
1424	tempfunc_result(:) = MAX(MIN( 1._r_std, tempfunc_result(:) ), zero)
1425
1426	END FUNCTION control_temp_func
1427
1428	!!!qcj++ peatland new functions
1429	!================================================================================================================================
1430
1431	FUNCTION control_moist_func_peat (npts, mc_peat) RESULT (moistfunc_result)
1432
1433	!! 0. Variable and parameter declaration
1434
1435	!! 0.1 Input variables
1436
1437	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixel (unitless)
1438	REAL(r_std), DIMENSION(npts), INTENT(in) :: mc_peat !! relative humidity (unitless)
1439	REAL(r_std),DIMENSION(40) :: mc !!!! used for Moyano et al., 2012, volumetric moisture, 0.01 interval
1440	REAL(r_std),DIMENSION(40) :: pcsr
1441	REAL(r_std),DIMENSION(40) :: sr
1442	REAL(r_std),DIMENSION(40) :: corgmat
1443	INTEGER(i_std) :: ind
1444	INTEGER(i_std) :: mc_ind
1445
1446	!! 0.2 Output variables
1447
1448	REAL(r_std), DIMENSION(npts) :: moistfunc_result !! Moisture control factor (0.25-1, unitless)
1449
1450	!! 0.3 Modified variables
1451	INTEGER :: ii,jj
1452	!! 0.4 Local variables
1453
1454	!_
1455	!================================================================================================================================
1456	DO jj=1,40
1457	mc(jj)=0.01+0.02*(jj-1)
1458	ENDDO
1459
1460	pcsr(:)=0.97509-0.48212mc(:)+1.83997(mc(:)*2)-1.56379(mc(:)**3)+ &
1461	0.098671.2+1.399440.05+0.179380.3-0.30307mc(:)1.2-0.30885mc(:)*0.3
1462
1463	DO jj=1,40
1464	IF (jj==1) THEN
1465	sr(jj) = pcsr(jj)
1466	ELSE
1467	sr(jj)=sr(jj-1)* pcsr(jj)
1468	ENDIF
1469	ENDDO
1470	sr(:)=sr(:)/ MAXVAL(sr)
1471
1472	corgmat(:)=sr(:)
1473	ind= MAXLOC(corgmat,1)
1474	corgmat(1:ind)=corgmat(1:ind)-MINVAL(corgmat(1:ind))
1475	corgmat(1:ind)=corgmat(1:ind)/MAXVAL(corgmat(1:ind))
1476
1477	DO ii=1,npts
1478	mc_ind = MIN(40, MAX(1, INT(mc_peat(ii)/0.02)+1))
1479	moistfunc_result(ii)= corgmat(mc_ind)
1480	moistfunc_result(ii)= MIN(un,MAX(zero,moistfunc_result(ii)))
1481	ENDDO
1482
1483	END FUNCTION control_moist_func_peat
1484	!_
1485	!================================================================================================================================
1486	!!!! NOTE: Here only include CASE(1) function for peat
1487	FUNCTION control_temp_func_peat (npts, temp_in, frozen_respiration_func) RESULT (tempfunc_result)
1488	!! 0. Variable and parameter declaration
1489
1490	!! 0.1 Input variables
1491	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of land pixels (unitless)
1492	REAL(r_std), DIMENSION(npts), INTENT(in) :: temp_in !! Temperature (K)
1493	INTEGER(i_std), INTENT(in) :: frozen_respiration_func
1494	!! 0.2 Output variables
1495	REAL(r_std), DIMENSION(npts) :: tempfunc_result !! Temperature control factor (0-1, unitless)
1496
1497	!! 0.3 Modified variables
1498
1499	!! 0.4 Local variables
1500
1501	!_
1502	!================================================================================================================================
1503
1504	WHERE (temp_in(:) .GT. ZeroCelsius ) !!! normal as above
1505	tempfunc_result(:) = EXP( log(q10_peat) * ( temp_in(:) - (ZeroCelsius+30.) ) / 10. )
1506	ELSEWHERE (temp_in(:) .GT. ZeroCelsius - 1. ) !!! linear dropoff to zero
1507	tempfunc_result(:) = (temp_in(:) - (ZeroCelsius - 1.)) * EXP( log(q10_peat) * ( ZeroCelsius - (ZeroCelsius+30.) ) / 10. )
1508	ELSEWHERE !!! zero
1509	tempfunc_result(:) = zero
1510	endwhere
1511
1512	tempfunc_result(:) = MAX(MIN( 1._r_std, tempfunc_result(:) ), zero)
1513
1514	END FUNCTION control_temp_func_peat
1515	!================================================================================================================================
1516
1517	END MODULE stomate_litter

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

Download in other formats: