Context Navigation

stomate_litter.f90 @ 7481

Last change on this file since 7481 was 6939, checked in by jinfeng.chang, 4 years ago
Update ORCHIDEE-GMv3.2 for publication
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 67.5 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_stomate/stomate_litter.f90 @ 7481

Download in other formats: