Context Navigation

stomate_litter.f90 @ 7346

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

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

Context Navigation

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

Download in other formats: