Context Navigation

source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_stomate/stomate_season.f90 @ 5816

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

Line
1	! =================================================================================================================================
2	! MODULE : stomate_season
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006). This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
7	!
8	!>\BRIEF This module calculates long-term meteorological parameters from daily temperatures
9	!! and precipitations (essentially for phenology).
10	!!
11	!!\n DESCRIPTION: None
12	!!
13	!! RECENT CHANGE(S): None
14	!!
15	!! SVN :
16	!! $HeadURL$
17	!! $Date$
18	!! $Revision$
19	!! \n
20	!_ =================================================================================================================================
21
22	MODULE stomate_season
23
24	! modules used:
25	USE xios_orchidee
26	USE ioipsl_para
27	USE stomate_data
28	USE constantes
29	USE constantes_soil
30	USE pft_parameters
31	USE solar
32	USE grid
33
34	IMPLICIT NONE
35
36	! private & public routines
37
38	PRIVATE
39	PUBLIC season,season_clear
40
41	LOGICAL, SAVE :: firstcall_season = .TRUE. !! first call (true/false)
42	!$OMP THREADPRIVATE(firstcall_season)
43
44	CONTAINS
45
46	!! ================================================================================================================================
47	!! SUBROUTINE : season_clear
48	!!
49	!>\BRIEF Flag setting
50	!!
51	!! DESCRIPTION : This subroutine sets flags ::firstcall_season, to .TRUE., and therefore activates
52	!! section 1.1 of the ::season subroutine which writes messages to the output. \n
53	!! This subroutine is called at the end of the subroutine ::stomate_clear, in the
54	!! module ::stomate.
55	!!
56	!! RECENT CHANGE(S):None
57	!!
58	!! MAIN OUTPUT VARIABLE(S): ::firstcall_season
59	!!
60	!! REFERENCE(S) : None
61	!!
62	!! FLOWCHART : None
63	!! \n
64	!_ =================================================================================================================================
65
66	SUBROUTINE season_clear
67	firstcall_season=.TRUE.
68	END SUBROUTINE season_clear
69
70
71
72	!! ================================================================================================================================
73	!! SUBROUTINE : season
74	!!
75	!>\BRIEF This subroutine calculates many of the long-term biometeorological variables
76	!! needed in the phenology subroutines and in the calculation of long-term vegetation
77	!! dynamics in the LPJ DGVM.
78	!!
79	!! DESCRIPTION This subroutine is called by the module ::stomate before LPJ, and mainly deals
80	!! with the calculation of long-term meteorological variables, carbon fluxes and
81	!! vegetation-related variables that are used to calculate vegetation dynamics
82	!! in the stomate modules relating to the phenology and to the longer-term
83	!! changes in vegetation type and fractional cover in the LPJ DGVM modules. \n
84	!! In sections 2 to 5, longer-term meteorological variables are calculated. The
85	!! long-term moisture availabilities are used in the leaf onset and senescence
86	!! phenological models ::stomate_phenology and ::stomate_turnover that require
87	!! a moisture condition. The long term temperatures are also required for phenology
88	!! but in addition they are used in calculations of C flux and the presence and
89	!! establishment of vegetation patterns on a longer timescale in the LPJ DGVM modules.
90	!! Finally the monthly soil humidity/relative soil moisture is used in the C
91	!! allocation module ::stomate_alloc. \n
92	!! Sections 12 to 14 also calculate long-term variables of C fluxes, including NPP,
93	!! turnover and GPP. The first two are used to calculate long-term vegetation
94	!! dynamics and land cover change in the LPJ DVGM modules. The weekly GPP is used
95	!! to determine the dormancy onset and time-length, as described below. \n
96	!! The long-term variables described above are are used in certain vegetation
97	!! dynamics processes in order to maintain consistency with the basic hypotheses of the
98	!! parameterisations of LPJ, which operates on a one year time step (Krinner et al., 2005).
99	!! In order to reduce the computer memory requirements, short-term variables (e.g. daily
100	!! temperatures) are not stored for averaging over a longer period. Instead
101	!! the long-term variables (e.g. monthly temperature) are calculated at each time step
102	!! using a linear relaxation method, following the equation:
103	!! \latexonly
104	!! \input{season_lin_relax_eqn1.tex}
105	!! \endlatexonly
106	!! \n
107	!! The long-term variables are therefore updated on a daily basis. This method allows for
108	!! smooth temporal variation of the long-term variables which are used to calculate
109	!! vegetation dynamics (Krinner et al., 2005). \n
110	!! Sections 6 to 11 calculate the variables required for to determine leaf onset in
111	!! the module ::stomate_phenology.
112	!! These include :
113	!! - the dormance onset and time-length, when GPP is below a certain threshold \n
114	!! - the growing degree days (GDD), which is the sum of daily temperatures greater than
115	!! -5 degrees C, either since the onset of the dormancy period, or since midwinter \n
116	!! - the number of chilling days, which is the number of days with a daily temperature
117	!! lower than a PFT-dependent threshold since the beginning of the dormancy period \n
118	!! - the number of growing days, which is the the number of days with a temperature
119	!! greater than -5 degrees C since the onset of the dormancy period \n
120	!! - the time since the minimum moisture availability in the dormancy period. \n
121	!! These variables are used to determine the conditions needed for the start of the
122	!! leaf growing season. The specific models to which they correspond are given below. \n
123	!! Sections 15 to 20 are used to update the maximum/minimum or the sum of various
124	!! meteorological and biological variables during the year that are required for
125	!! calculating either leaf onset or longer-term vegetation dynamics the following year. \n
126	!! At the end of the year, these variables are updated from "thisyear" to "lastyear",
127	!! in Section 21 of this subroutine, for use the following year. \n
128	!! Finally the probably amount of herbivore consumption is calculated in Section 22,
129	!! following McNaughton et al. (1989).
130	!!
131	!! RECENT CHANGE(S): None
132	!!
133	!! MAIN OUTPUT VARIABLE(S): :: herbivores,
134	!! :: maxmoiavail_lastyear, :: maxmoiavail_thisyear,
135	!! :: minmoiavail_lastyear, :: minmoiavail_thisyear,
136	!! :: maxgppweek_lastyear, :: maxgppweek_thisyear,
137	!! :: gdd0_lastyear, :: gdd0_thisyear,
138	!! :: precip_lastyear, :: precip_thisyear,
139	!! :: lm_lastyearmax, :: lm_thisyearmax,
140	!! :: maxfpc_lastyear, :: maxfpc_thisyear,
141	!! :: moiavail_month, :: moiavail_week,
142	!! :: t2m_longterm, :: t2m_month, :: t2m_week,
143	!! :: tsoil_month, :: soilhum_month,
144	!! :: npp_longterm, :: turnover_longterm, :: gpp_week,
145	!! :: gdd_m5_dormance, :: gdd_midwinter,
146	!! :: ncd_dormance, :: ngd_minus5, :: time_lowgpp,
147	!! :: time_hum_min, :: hum_min_dormance
148	!!
149	!! REFERENCES :
150	!! - Krinner, G., N. Viovy, N. de Noblet-Ducoudre, J. Ogee, J. Polcher, P.
151	!! Friedlingstein, P. Ciais, S. Sitch and I.C. Prentice (2005), A dynamic global
152	!! vegetation model for studies of the coupled atmosphere-biosphere system, Global
153	!! Biogeochemical Cycles, 19, doi:10.1029/2003GB002199.
154	!! - McNaughton, S.J., M. Oesterheld, D.A. Frank and K.J. Williams (1989),
155	!! Ecosystem-level patterns of primary productivity and herbivory in terrestrial
156	!! habitats, Nature, 341, 142-144.
157	!!
158	!! FLOWCHART :
159	!! \latexonly
160	!! \includegraphics[scale = 1]{season_flowchart_part1.png}
161	!! \includegraphics[scale = 1]{season_flowchart_part2.png}
162	!! \includegraphics[scale = 1]{season_flowchart_part3.png}
163	!! \endlatexonly
164	!! \n
165	!_ =================================================================================================================================
166
167	SUBROUTINE season (npts, dt, EndOfYear, veget, veget_max, &
168	moiavail_daily, t2m_daily, tsoil_daily, soilhum_daily, lalo, &
169	precip_daily, npp_daily, biomass, turnover_daily, gpp_daily, when_growthinit, &
170	maxmoiavail_lastyear, maxmoiavail_thisyear, &
171	minmoiavail_lastyear, minmoiavail_thisyear, &
172	maxgppweek_lastyear, maxgppweek_thisyear, &
173	gdd0_lastyear, gdd0_thisyear, &
174	precip_lastyear, precip_thisyear, &
175	lm_lastyearmax, lm_thisyearmax, &
176	maxfpc_lastyear, maxfpc_thisyear, &
177	moiavail_month, moiavail_week, t2m_longterm, tau_longterm, t2m_month, t2m_week, &
178	tsoil_month, soilhum_month, &
179	npp_longterm, turnover_longterm, gpp_week, &
180	gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
181	time_hum_min, hum_min_dormance, gdd_init_date , &
182	gdd_from_growthinit, herbivores, &
183	Tseason, Tseason_length, Tseason_tmp, &
184	Tmin_spring_time, t2m_min_daily, begin_leaves, onset_date, &
185	!gmjc
186	t2m_14, sla_calc)
187	!end gmjc
188
189	!
190	!! 0. Variable and parameter declaration
191	!
192
193	!
194	!! 0.1 Input variables
195	!
196	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid cells (unitless)
197	REAL(r_std), INTENT(in) :: dt !! time step in days (dt_days)
198	LOGICAL, INTENT(in) :: EndOfYear !! update yearly variables? (true/false)
199	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget !! coverage fraction of a PFT. Here: fraction of
200	!! total ground. (0-1, unitless)
201	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT (for LAI ->
202	!! infinity) (0-1, unitless)Here: fraction of
203	!! total ground.
204	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: moiavail_daily !! Daily moisture availability (0-1, unitless)
205	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_daily !! Daily 2 meter temperature (K)
206	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: tsoil_daily !! Daily soil temperature (K)
207	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: soilhum_daily !! Daily soil humidity (0-1, unitless)
208	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: lalo !! array of lat/lon
209	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_daily !! Daily mean precipitation @tex ($mm day^{-1}$)
210	!! @endtex
211	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: npp_daily !! daily net primary productivity @tex ($gC m^{-2}
212	!! day^{-1}$) @endtex
213	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2} of ground$) @endtex
214	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: turnover_daily !! Turnover rates @tex ($gC m^{-2} day^{-1}$)
215	!! @endtex
216	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: gpp_daily !! daily gross primary productivity
217	!! (Here: @tex $gC m^{-2} of total ground
218	!! day^{-1}$) @endtex
219	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: when_growthinit !! how many days ago was the beginning of the
220	!! growing season (days)
221	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: begin_leaves !! signal to start putting leaves on (true/false)
222	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_min_daily !! Daily minimum 2-meter temperature (K)
223
224	!
225	!! 0.2 Output variables
226	! (diagnostic)
227	!
228	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: herbivores !! time constant of probability of a leaf to be
229	!! eaten by a herbivore (days)
230
231	!
232	!! 0.3 Modified variables
233	!
234	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxmoiavail_lastyear !! last year's maximum moisture
235	!! availability (0-1, unitless)
236	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxmoiavail_thisyear !! this year's maximum moisture
237	!! availability (0-1, unitless)
238	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: minmoiavail_lastyear !! last year's minimum moisture
239	!! availability (0-1, unitless)
240	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: minmoiavail_thisyear !! this year's minimum moisture
241	!! availability (0-1, unitless)
242	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxgppweek_lastyear !! last year's maximum weekly GPP
243	!! @tex ($gC m^{-2} week^{-1}$) @endtex
244	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxgppweek_thisyear !! this year's maximum weekly GPP
245	!! @tex ($gC m^{-2} week^{-1}$) @endtex
246	REAL(r_std), DIMENSION(npts), INTENT(inout) :: gdd0_lastyear !! last year's annual GDD0 (C)
247	REAL(r_std), DIMENSION(npts), INTENT(inout) :: gdd0_thisyear !! this year's annual GDD0 (C)
248	REAL(r_std), DIMENSION(npts), INTENT(inout) :: precip_lastyear !! last year's annual precipitation
249	!! @tex ($mm year^{-1}$) @endtex
250	REAL(r_std), DIMENSION(npts), INTENT(inout) :: precip_thisyear !! this year's annual precipitation
251	!! @tex ($mm year^{-1}$) @endtex
252	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass, for each
253	!! PFT @tex ($gC m^{-2}$) @endtex
254	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_thisyearmax !! this year's maximum leaf mass, for each
255	!! PFT @tex ($gC m^{-2}$) @endtex
256	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxfpc_lastyear !! last year's maximum fpc for each PFT, on
257	!! ground (0-1, unitless)
258	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxfpc_thisyear !! this year's maximum fpc for each PFT, on
259	!! ground (0-1, unitless)
260	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_month !! "monthly" moisture availability
261	!! (0-1, unitless)
262	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_week !! "weekly" moisture availability
263	!! (0-1, unitless)
264	REAL(r_std), DIMENSION(npts), INTENT(inout) :: t2m_longterm !! "long term" 2-meter temperatures (K)
265	REAL(r_std), INTENT(inout) :: tau_longterm
266	REAL(r_std), DIMENSION(npts), INTENT(inout) :: t2m_month !! "monthly" 2-meter temperatures (K)
267	REAL(r_std), DIMENSION(npts), INTENT(inout) :: t2m_week !! "weekly" 2-meter temperatures (K)
268	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Tseason !! "seasonal" 2-meter temperatures (K),
269	!! used to constrain boreal treeline
270	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Tseason_tmp !! temporary variable to calculate Tseason
271	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Tseason_length !! temporary variable to calculate Tseason:
272	!! number of days when t2m_week is higher than 0 degree
273	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: Tmin_spring_time !! Number of days after begin_leaves (leaf onset)
274	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: onset_date !! Date in the year at when the leaves started to grow(begin_leaves)
275	REAL(r_std), DIMENSION(npts,nbdl), INTENT(inout) :: tsoil_month !! "monthly" soil temperatures (K)
276	REAL(r_std), DIMENSION(npts,nbdl), INTENT(inout) :: soilhum_month !! "monthly" soil humidity (0-1, unitless)
277	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: npp_longterm !! "long term" net primary productivity
278	!! @tex ($gC m^{-2} year^{-1}$) @endtex
279	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: turnover_longterm !! "long term" turnover rate
280	!! @tex ($gC m^{-2} year^{-1}$) @endtex
281	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gpp_week !! "weekly" GPP @tex ($gC m^{-2} day^{-1}$)
282	!! @endtex
283	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_m5_dormance !! growing degree days above threshold -5
284	!! deg. C (C)
285	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_midwinter !! growing degree days since midwinter (C)
286	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ncd_dormance !! number of chilling days since leaves
287	!! were lost (days)
288	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ngd_minus5 !! number of growing days (days)
289	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: time_hum_min !! time elapsed since strongest moisture
290	!! availability (days)
291	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: hum_min_dormance !! minimum moisture during dormance
292	!! (0-1, unitless)
293	REAL(r_std), DIMENSION(npts,2), INTENT(inout) :: gdd_init_date !! inital date for gdd count
294	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_from_growthinit !! growing degree days, threshold 0 deg. C
295	!! since beginning of season
296	!gmjc
297	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: sla_calc
298	! "14days" 2-meter temperatures (K)
299	REAL(r_std), DIMENSION(npts), INTENT(inout) :: t2m_14
300	!end gmjc
301
302	!
303	!! 0.4 Local variables
304	!
305	INTEGER(i_std) :: j !! indices (unitless)
306	REAL(r_std) :: ncd_max !! maximum ncd (to avoid floating point
307	!! underflows) (days)
308	REAL(r_std), DIMENSION(npts) :: sumfpc_nat !! sum of natural fpcs (0-1, unitless)
309	!! [DISPENSABLE]
310	REAL(r_std), DIMENSION(npts,nvm) :: weighttot !! weight of biomass @tex ($gC m^{-2}$)
311	!! @endtex
312	REAL(r_std), DIMENSION(npts,nvm) :: nlflong_nat !! natural long-term leaf NPP
313	!! @tex ($gC m^{-2} year^{-1}$) @endtex
314	REAL(r_std), DIMENSION(npts,nvm) :: green_age !! residence time of green tissue (years)
315	REAL(r_std), DIMENSION(npts) :: consumption !! herbivore consumption
316	!! @tex ($gC m^{-2} day^{-1}$) @endtex
317	REAL(r_std), DIMENSION(npts) :: fracnat !! fraction of each gridcell occupied by
318	!! natural vegetation (0-1, unitless)
319	REAL(r_std), DIMENSION(npts) :: solad !! maximal radation during current day
320	!! (clear sky condition)
321	REAL(r_std), DIMENSION(npts) :: solai !! maximal radation during current day
322	!! (clear sky condition)
323	REAL(r_std), DIMENSION(npts) :: cloud !! cloud fraction
324
325	!_ =================================================================================================================================
326
327	IF (printlev>=3) WRITE(numout,*) 'Entering season'
328
329	!
330	!! 1. Initializations
331	!! 1.1 Calculate ::ncd_max - the maximum possible NCD (number of chilling days) as:
332	!! \latexonly
333	!! \input{season_ncdmax_eqn2.tex}
334	!! \endlatexonly
335	!! \n
336	!! where one_year is 1 year in seconds (defined in ::constantes).
337	!
338	ncd_max = ncd_max_year * one_year
339
340	IF ( firstcall_season ) THEN
341
342	!
343	!! 1.2 first call - output message giving the setting of the ::gppfrac_dormance,
344	!! ::hvc1, ::hvc2 and ::leaf_frac (as a percentage) parameters which are
345	!! defined at the beginning of this subroutine. Also outputs the value of
346	!! ::ncd_max.
347	!
348
349	IF ( printlev>=1 ) THEN
350
351	WRITE(numout,*) 'season: '
352
353	WRITE(numout,) ' > maximum GPP/GGP_max ratio for dormance (::gppfrac_dormance) :',gppfrac_dormance !
354
355	WRITE(numout,) ' > maximum possible ncd(days) (::ncd_max) : ',ncd_max !
356
357	WRITE(numout,) ' > herbivore consumption C (gC/m2/d) as a function of NPP (gC/m2/d): (::hvc1) (::hvc2)' !
358	WRITE(numout,) ' C=',hvc1,' NPP^',hvc2
359	WRITE(numout,) ' > for herbivores, suppose that (::leaf_frac_hvc)',leaf_frac_hvc100., &
360	'% of NPP is allocated to leaves'
361
362	ENDIF
363
364	!
365	!! 1.3 Check whether longer-term meteorological and phenology-related variables are not initialized
366	!! to less than zero. If the following meteorological variables are less than ::min_stomate which is set
367	!! to 1.E-8 in ::grid, then they are set to the current daily value. If the phenology-
368	!! related variables are less than ::min_stomate, they are set to "undef".
369	!! Warning messages are output in each case, and are also output if "long-term" C fluxes GPP,
370	!! NPP and turnover are less than ::min_stomate.
371	!
372
373	!! 1.3.1 moisture availabilities (i.e. relative soil moisture, including the root soil profile):
374
375	!! 1.3.1.1 "monthly" (::moiavail_month)
376	!MM PAS PARALLELISE!!
377	IF ( ABS( SUM( moiavail_month(:,2:nvm) ) ) .LT. min_stomate ) THEN
378
379	! in this case, set them it today's moisture availability
380	WRITE(numout,*) 'Warning! We have to initialize the ''monthly'' moisture availabilities. '
381	moiavail_month(:,:) = moiavail_daily(:,:)
382
383	ENDIF
384
385	!! 1.3.1.2 "weekly" (::moiavail_week)
386
387	IF ( ABS( SUM( moiavail_week(:,2:nvm) ) ) .LT. min_stomate ) THEN
388
389	! in this case, set them it today's moisture availability
390	WRITE(numout,*) 'Warning! We have to initialize the ''weekly'' moisture availabilities. '
391	moiavail_week(:,:) = moiavail_daily(:,:)
392
393	ENDIF
394
395	!! 1.3.2 2-meter temperatures:
396
397	!! "monthly" (::t2m_month)
398	IF ( ABS( SUM( t2m_month(:) ) ) .LT. min_stomate ) THEN
399
400	! in this case, set them to today's temperature
401	WRITE(numout,*) 'Warning! We have to initialize the ''monthly'' 2m temperatures.'
402	t2m_month(:) = t2m_daily(:)
403
404	ENDIF
405
406	!! "weekly" (::2m_week)
407	IF ( ABS( SUM( t2m_week(:) ) ) .LT. min_stomate ) THEN
408
409	! in this case, set them to today's temperature
410	WRITE(numout,*) 'Warning! We have to initialize the ''weekly'' 2m temperatures.'
411	t2m_week(:) = t2m_daily(:)
412
413	ENDIF
414
415	!gmjc
416	!! 1.3.2.2 "14days" (::t2m_14)
417
418	IF ( ABS( SUM( t2m_14(:) ) ) .LT. min_stomate ) THEN
419
420	! in this case, set them to today's temperature
421	WRITE(numout,*) 'Warning! We have to initialize the ''14days'' 2m temperatures.'
422	t2m_14(:) = t2m_daily(:)
423
424	ENDIF
425	!end gmjc
426
427	!! 1.3.3 "monthly" soil temperatures (::tsoil_month):
428	!MM PAS PARALLELISE!!
429	IF ( ABS( SUM( tsoil_month(:,:) ) ) .LT. min_stomate ) THEN
430
431	! in this case, set them to today's temperature
432	WRITE(numout,*) 'Warning!'// &
433	' We have to initialize the ''monthly'' soil temperatures.'
434	tsoil_month(:,:) = tsoil_daily(:,:)
435
436	ENDIF
437
438	!! 1.3.4 "monthly" soil humidity (::soilhum_month):
439	IF ( ABS( SUM( soilhum_month(:,:) ) ) .LT. min_stomate ) THEN
440
441	! in this case, set them to today's humidity
442	WRITE(numout,*) 'Warning!'// &
443	' We have to initialize the ''monthly'' soil humidity.'
444	soilhum_month(:,:) = soilhum_daily(:,:)
445
446	ENDIF
447
448	!! 1.3.5 growing degree days, threshold -5 deg C (::gdd_m5_dormance):
449	IF ( ABS( SUM( gdd_m5_dormance(:,2:nvm) ) ) .LT. min_stomate ) THEN
450	WRITE(numout,*) 'Warning! Growing degree days (-5 deg) are initialized to ''undef''.'
451	gdd_m5_dormance(:,:) = undef
452	ENDIF
453
454	!! 1.3.6 growing degree days since midwinter (::gdd_midwinter):
455	IF ( ABS( SUM( gdd_midwinter(:,2:nvm) ) ) .LT. min_stomate ) THEN
456	WRITE(numout,*) 'Warning! Growing degree days since midwinter' // &
457	' are initialized to ''undef''.'
458	gdd_midwinter(:,:) = undef
459	ENDIF
460
461	!! 1.3.7 number of chilling days since leaves were lost (::ncd_dormance):
462	IF ( ABS( SUM( ncd_dormance(:,2:nvm) ) ) .LT. min_stomate ) THEN
463	WRITE(numout,*) 'Warning! Number of chilling days is initialized to ''undef''.'
464	ncd_dormance(:,:) = undef
465	ENDIF
466
467	!! 1.3.8 number of growing days, threshold -5 deg C (::ngd_minus5):
468	IF ( ABS( SUM( ngd_minus5(:,2:nvm) ) ) .LT. min_stomate ) THEN
469	WRITE(numout,*) 'Warning! Number of growing days (-5 deg) is initialized to 0.'
470	ENDIF
471
472	!! 1.3.9 "long term" npp (::npp_longterm):
473	IF ( ABS( SUM( npp_longterm(:,2:nvm) ) ) .LT. min_stomate ) THEN
474	WRITE(numout,*) 'Warning! Long term NPP is initialized to 0.'
475	ENDIF
476
477	!! 1.3.10 "long term" turnover (::turnover_longterm):
478	IF ( ABS( SUM( turnover_longterm(:,2:nvm,:,:) ) ) .LT. min_stomate ) THEN
479	WRITE(numout,*) 'Warning! Long term turnover is initialized to 0.'
480	ENDIF
481
482	!! 1.3.11 "weekly" GPP (::gpp_week):
483	IF ( ABS( SUM( gpp_week(:,2:nvm) ) ) .LT. min_stomate ) THEN
484	WRITE(numout,*) 'Warning! Weekly GPP is initialized to 0.'
485	ENDIF
486
487	!! 1.3.12 minimum moisture availabilities (::minmoiavail_thisyear)
488	!! In this case, if the minimum moisture
489	!! is less than ::min_stomate, the variable is set to ::large_value,
490	!! which is defined in ::stomate_constants as 1.E33_r_std:
491	IF ( ABS( SUM( minmoiavail_thisyear(:,2:nvm) ) ) .LT. min_stomate ) THEN
492
493	! in this case, set them to a very high value
494	WRITE(numout,*) 'Warning! We have to initialize this year''s minimum '// &
495	'moisture availabilities.'
496	minmoiavail_thisyear(:,:) = large_value
497
498	ENDIF
499
500	!
501	!! 1.4 reset firstcall_season flag
502	!
503
504	firstcall_season = .FALSE.
505
506	ENDIF
507
508	! determine min yearly clear sky radiation (solstice) as beginning for gdd count
509	cloud(:) = zero
510	CALL downward_solar_flux (npts, lalo(:,1),calendar_str,julian_diff,12.,cloud,1,solad,solai)
511
512	WHERE (solad(:) .LT. gdd_init_date(:,2))
513	gdd_init_date(:,1)= julian_diff
514	gdd_init_date(:,2)= solad(:)
515	ENDWHERE
516
517
518	!
519	!! NOTE: Sections 2. to 5. compute slowly-varying, "long-term" (i.e. weekly/monthly)
520	!! input variables using a linear relaxation method, following the equation:
521	!! \latexonly
522	!! \input{season_lin_relax_eqn1.tex}
523	!! \endlatexonly
524	!! \n
525	!! as described in the introduction to this subroutine (see Krinner et al., 2005).
526	!! The time constant in the above equation is given for each of the variables
527	!! described.
528	!
529
530	!
531	!! 2. Moisture availability (relative soil moisture, including the root profile) :
532	!! The time constants (as in the above equation) for these calculations are
533	!! given by the parameters ::tau_hum_month (for the monthly
534	!! calculation) or ::tau_hum_week (for the weekly),
535	!! which are set in ::stomate_data to be 20 and 7 days respectively.
536	!! If the moisture availability is less than zero, it is set to zero.
537	!! These variables are mostly used in the phenological leaf onset and senescence
538	!! models in the modules ::stomate_phenology and ::stomate_turnover,
539	!! respectively. They are used in models which require a moisture limitation
540	!! condition. These are the 'hum', 'moi', 'humgdd' and 'moigdd' onset models,
541	!! and the 'mixed' and 'dry' (weekly only) senescence models. In addition,
542	!! the weekly moisture availability is used to calculate the limitation
543	!! on the fraction of NPP allocated to different compartments in the module
544	!! ::stomate_alloc.
545	!
546
547	!
548	!! 2.1 "monthly" (::moiavail_month)
549	!
550
551	moiavail_month = ( moiavail_month * ( tau_hum_month - dt ) + &
552	moiavail_daily * dt ) / tau_hum_month
553
554	DO j = 2,nvm ! Loop over # PFTs
555	WHERE ( ABS(moiavail_month(:,j)) .LT. EPSILON(zero) )
556	moiavail_month(:,j) = zero
557	ENDWHERE
558	ENDDO
559
560	!
561	!! 2.2 "weekly" (::moiavail_week)
562	!
563
564	moiavail_week = ( moiavail_week * ( tau_hum_week - dt ) + &
565	moiavail_daily * dt ) / tau_hum_week
566
567	DO j = 2,nvm ! Loop over # PFTs
568	WHERE ( ABS(moiavail_week(:,j)) .LT. EPSILON(zero) )
569	moiavail_week(:,j) = zero
570	ENDWHERE
571	ENDDO
572
573	!
574	!! 3. 2-meter temperatures
575	!! The time constants for the "long-term", "monthly" and "weekly" 2-meter
576	!! temperatures are given by the parameters ::tau_longterm_max,
577	!! ::tau_t2m_month, and ::tau_t2m_week,
578	!! which are set in ::stomate_data to be 3 * one year (in seconds, as
579	!! described above) and 20 and 7 days respectively.
580	!! If the temperature is less than zero, it is set to zero.
581	!! In addition the "long term" temperature is written to the history file. \n
582	!! These variables are used in many different modules of STOMATE.
583	!! The longterm t2m is limied to the intervall [tlong_ref_min, tlong_ref_max]
584	!! using the equation:
585	!! \latexonly
586	!! \input{season_t2m_long_ref_eqn3.tex}
587	!! \endlatexonly
588	!! \n
589	!! The "monthly" and "weekly" temperature is used in the modules
590	!! ::stomate_phenology and ::stomate_turnover for the onset and senescence
591	!! phenological models which require a temperature condition. In addition
592	!! the "monthly" temperature is used in ::lpj_constraints to determine
593	!! the presence and regeneration of the vegetation and in the module
594	!! ::stomate_assimtemp to calculate photosynthesis temperatures.
595
596	! Update tau_longterm
597	tau_longterm = MIN(tau_longterm+dt,tau_longterm_max)
598
599	! Recalculate the reference temperature using the old reference temperature and the current temperature
600	t2m_longterm(:) = ( t2m_longterm(:) * ( tau_longterm - dt ) + &
601	t2m_daily(:) * dt ) / tau_longterm
602
603	! The longterm reference is not allowed to go outside the interval [tlong_ref_min, tlong_ref_max]
604	t2m_longterm(:) = MAX( tlong_ref_min, MIN( tlong_ref_max, t2m_longterm(:) ) )
605
606	CALL xios_orchidee_send_field("T2M_LONGTERM",t2m_longterm)
607
608	CALL histwrite_p (hist_id_stomate, 'T2M_LONGTERM', itime, &
609	t2m_longterm, npts, hori_index)
610
611	!
612	!! 3.3 "monthly" (::t2m_month)
613	!
614
615	t2m_month = ( t2m_month * ( tau_t2m_month - dt ) + &
616	t2m_daily * dt ) / tau_t2m_month
617
618	WHERE ( ABS(t2m_month(:)) .LT. EPSILON(zero) )
619	t2m_month(:) = zero
620	ENDWHERE
621
622	! Store the day in onset_date when the leaves started grow. julian_diff is the day in the year [1-366].
623	WHERE ( begin_leaves(:,:) )
624	onset_date(:,:)=julian_diff
625	ELSEWHERE
626	onset_date(:,:)=0
627	ENDWHERE
628
629	! Calculate "seasonal" temperature
630	WHERE ( t2m_week (:) .GT. ZeroCelsius )
631	Tseason_tmp(:) = Tseason_tmp(:) + t2m_week(:)
632	Tseason_length(:)=Tseason_length(:) + dt
633	ENDWHERE
634
635	! calculate Tmin in spring (after onset)
636	DO j=1, nvm
637	IF (leaf_tab(j)==1 .AND. pheno_type(j)==2) THEN
638	! leaf_tab=broadleaf and pheno_typ=summergreen
639	! Only treat broadleag and summergreen pfts
640
641	WHERE ( Tmin_spring_time(:,j)>0 .AND. (Tmin_spring_time(:,j)<spring_days_max) )
642	Tmin_spring_time(:,j)=Tmin_spring_time(:,j)+1
643	ELSEWHERE
644	Tmin_spring_time(:,j)=0
645	ENDWHERE
646
647	WHERE ( begin_leaves(:,j) )
648	Tmin_spring_time(:,j)=1
649	ENDWHERE
650	END IF
651	END DO
652
653	!
654	!! 3.4 "weekly" (::t2m_week)
655	!
656
657	t2m_week = ( t2m_week * ( tau_t2m_week - dt ) + &
658	t2m_daily * dt ) / tau_t2m_week
659
660	WHERE ( ABS(t2m_week(:)) .LT. EPSILON(zero) )
661	t2m_week(:) = zero
662	ENDWHERE
663	!gmjc
664	!
665	! 3.5 "14days"
666	!
667
668	t2m_14 = ( t2m_14 * ( tau_t2m_14 - dt ) + &
669	t2m_daily * dt ) / tau_t2m_14
670
671	WHERE ( ABS(t2m_14(:)) .LT. EPSILON(zero) )
672	t2m_14(:) = t2m_daily(:)
673	ENDWHERE
674	!end gmjc
675	!
676	!! 4. "monthly" soil temperatures (::tsoil_month)
677	!! The time constant is given by the parameter ::tau_tsoil_month,
678	!! which is set in ::stomate_data to be 20 days.
679	!! If the monthly soil temperature is less than zero, it is set to zero. \n
680	!! This variable is used in the ::stomate_allocation module.
681	!
682
683	tsoil_month = ( tsoil_month * ( tau_tsoil_month - dt ) + &
684	tsoil_daily(:,:) * dt ) / tau_tsoil_month
685
686	WHERE ( ABS(tsoil_month(:,:)) .LT. EPSILON(zero) )
687	tsoil_month(:,:) = zero
688	ENDWHERE
689
690	!
691	!! 5. ''monthly'' soil humidity (or relative soil moisture - ::soilhum_month)
692	!! The time constant is given by the parameter ::tau_soilhum_month,
693	!! which is set in ::stomate_data to be 20 days.
694	!! If the monthly soil humidity is less than zero, it is set to zero.
695	!
696
697	soilhum_month = ( soilhum_month * ( tau_soilhum_month - dt ) + &
698	soilhum_daily * dt ) / tau_soilhum_month
699
700	WHERE ( ABS(soilhum_month(:,:)) .LT. EPSILON(zero) )
701	soilhum_month(:,:) = zero
702	ENDWHERE
703
704	!
705	!! 6. Calculate the dormance time-length (::time_lowgpp).
706	!! The dormancy time-length is increased by the stomate time step when
707	!! the weekly GPP is below one of two thresholds,
708	!! OR if the growing season started more than 2 years ago AND the amount of biomass
709	!! in the carbohydrate reserve is 4 times as much as the leaf biomass.
710	!! i.e. the plant is declared dormant if it is accumulating carbohydrates but not
711	!! using them, so that the beginning of the growing season can be detected .
712	! This last condition was added by Nicolas Viovy.
713	!! - NV: special case (3rd condition)
714	!! Otherwise, set ::time_lowgpp to zero. \n
715	!! The weekly GPP must be below either the parameter ::min_gpp_allowed, which is set at
716	!! the beginning of this subroutine to be 0.3gC/m^2/year,
717	!! OR it must be less than last year's maximum weekly GPP multiplied by the
718	!! maximal ratio of GPP to maximum GPP (::gppfrac_dormance), which is set to 0.2 at the
719	!! beginning of this subroutine. \n
720	!! This variable is used in the ::stomate_phenology module to allow the detection of the
721	!! beginning of the vegetation growing season for different onset models.
722	!! Each PFT is looped over.
723	!
724	!NVMODIF
725	!!$ DO j = 2,nvm ! Loop over # PFTs
726	!!$ WHERE ( ( gpp_week(:,j) .LT. min_gpp_allowed ) .OR. &
727	!!$ ( gpp_week(:,j) .LT. gppfrac_dormance * maxgppweek_lastyear(:,j) ) .OR. &
728	!!$ ( ( when_growthinit(:,j) .GT. 2.*one_year ) .AND. &
729	!!$ ( biomass(:,j,icarbres,icarbon) .GT. biomass(:,j,ileaf,icarbon)*4. ) ) )
730	!!$ ! WHERE ( ( gpp_week(:,j) .EQ. zero ) .OR. &
731	!!$ ! ( gpp_week(:,j) .LT. gppfrac_dormance * maxgppweek_lastyear(:,j) ) .OR. &
732	!!$ ! ( ( when_growthinit(:,j) .GT. 2.*one_year ) .AND. &
733	!!$ ! ( biomass(:,j,icarbres,icarbon) .GT. biomass(:,j,ileaf,icarbon)*4. ) ) )
734	!!$
735	!!$ time_lowgpp(:,j) = time_lowgpp(:,j) + dt
736	!!$
737	!!$ ELSEWHERE
738	!!$
739	!!$ time_lowgpp(:,j) = zero
740	!!$
741	!!$ ENDWHERE
742	!!$ ENDDO
743
744	!
745	!! 7. Calculate the growing degree days (GDD - ::gdd_m5_dormance).
746	!! This variable is the GDD sum of the temperatures higher than -5 degrees C.
747	!! It is inititalised to 0 at the beginning of the dormancy period (i.e.
748	!! when ::time_lowgpp>0), and is set to "undef" when there is GPP
749	!! (::time_lowgpp=0), as it is not used during the growing season. \n
750	!! ::gdd_m5_dormance is further scaled by (::tau_gdd -dt / ::tau_gdd),
751	!! - ::tau_gdd is set to 40. in the module ::stomate_constants. \n
752	! Nicolas Viovy - not sure why this is...
753	!! This variable is used in the ::stomate_phenology module for the
754	!! 'humgdd' and 'moigdd' leaf onset models. \n
755	!! Each PFT is looped over but ::gdd_m5_dormance is only calculated for
756	!! those PFTs for which a critical GDD is defined, i.e. those which are
757	!! assigned to the 'humgdd' or 'moigdd' models. \n
758	!! Finally if GDD sum is less than zero, then it is set to zero.
759	!
760
761	DO j = 2,nvm ! Loop over # PFTs
762
763	IF (.NOT. natural(j)) THEN
764	! reset counter: start of the growing season
765	WHERE ((when_growthinit(:,j) .EQ. zero))
766	gdd_from_growthinit(:,j) = zero
767	ENDWHERE
768	! increase gdd counter
769	WHERE ( t2m_daily(:) .GT. (ZeroCelsius) )
770	gdd_from_growthinit(:,j) = gdd_from_growthinit(:,j) + &
771	dt * ( t2m_daily(:) - (ZeroCelsius) )
772	ENDWHERE
773	ELSE
774	gdd_from_growthinit(:,j) = undef
775	ENDIF
776
777	! only for PFTs for which critical gdd is defined
778	! gdd_m5_dormance is set to 0 at the end of the growing season. It is set to undef
779	! at the beginning of the growing season.
780
781	IF ( ALL(pheno_gdd_crit(j,:) .NE. undef) ) THEN
782
783	!
784	!! 7.1 set to zero if undefined and there is no GPP
785	!
786
787	WHERE (gdd_init_date(:,1) .EQ. julian_diff)
788
789	gdd_m5_dormance(:,j) = zero
790
791	ENDWHERE
792
793	!
794	!! 7.2 set to undef if there is GPP
795	!
796
797	WHERE ( when_growthinit(:,j) .EQ. zero )
798
799	gdd_m5_dormance(:,j) = undef
800
801	ENDWHERE
802
803	!
804	!! 7.3 normal update as described above where ::gdd_m5_dormance is defined (not set to "undef").
805	!
806
807	WHERE ( ( t2m_daily(:) .GT. ( ZeroCelsius - gdd_threshold ) ) .AND. &
808	( gdd_m5_dormance(:,j) .NE. undef ) )
809	gdd_m5_dormance(:,j) = gdd_m5_dormance(:,j) + &
810	dt * ( t2m_daily(:) - ( ZeroCelsius - gdd_threshold ) )
811	ENDWHERE
812
813	WHERE ( gdd_m5_dormance(:,j) .NE. undef )
814	gdd_m5_dormance(:,j) = gdd_m5_dormance(:,j) * &
815	( tau_gdd - dt ) / tau_gdd
816	ENDWHERE
817
818	ENDIF
819
820	ENDDO
821
822	!
823	!! 7.4 Set to zero if GDD is less than zero.
824
825	DO j = 2,nvm ! Loop over # PFTs
826	WHERE ( ABS(gdd_m5_dormance(:,j)) .LT. EPSILON(zero) )
827	gdd_m5_dormance(:,j) = zero
828	ENDWHERE
829	ENDDO
830
831	!
832	!! 8. Calculate the growing degree days (GDD) since midwinter (::gdd_midwinter)
833	!! This variable represents the GDD sum of temperatures higher than a PFT-dependent
834	!! threshold (::ncdgdd_temp), since midwinter.
835	!! Midwinter is detected if the monthly temperature (::t2m_month) is lower than the weekly
836	!! temperature (::t2m_week) AND if the monthly temperature is lower than the long-term
837	!! temperature (t2m_longterm). These variables were calculated earlier in this subroutine. \n
838	!! ::gdd_midwinter is initialised to 0.0 when midwinter is detected, and from then on
839	!! increased with each temperature greater than ::ncdgdd_temp, which is defined
840	!! in the module ::stomate_constants. \n
841	!! ::gdd_midwinter is set to "undef" when midsummer is detected, following the opposite
842	!! conditions to those used to define midwinter. \n
843	!! The variable is used in the ::stomate_phenology module for the leaf onset model 'ncdgdd'.
844	!! Each PFT is looped over but the ::gdd_midwinter is only calculated for those
845	!! PFTs for which a critical 'ncdgdd' temperature is defined, i.e. those which are
846	!! assigned to the 'ncdgdd' model. \n
847	!
848
849	DO j = 2,nvm ! Loop over # PFTs
850
851	! only for PFTs for which ncdgdd_crittemp is defined
852
853	IF ( ncdgdd_temp(j) .NE. undef ) THEN
854
855	!
856	!! 8.1 set to 0 if undef and if we detect "midwinter"
857	!
858
859	!!$ WHERE ( ( gdd_midwinter(:,j) .EQ. undef ) .AND. &
860	!!$ ( t2m_month(:) .LT. t2m_week(:) ) .AND. &
861	!!$ ( t2m_month(:) .LT. t2m_longterm(:) ) )
862	WHERE (gdd_init_date(:,1) .EQ. julian_diff)
863
864	gdd_midwinter(:,j) = zero
865
866	ENDWHERE
867
868	!
869	!! 8.2 set to undef if we detect "midsummer"
870	!
871
872	WHERE ( ( t2m_month(:) .GT. t2m_week(:) ) .AND. &
873	( t2m_month(:) .GT. t2m_longterm(:) ) )
874
875	gdd_midwinter(:,j) = undef
876
877	ENDWHERE
878
879	!
880	!! 8.3 normal update as described above
881	!
882
883	WHERE ( ( gdd_midwinter(:,j) .NE. undef ) .AND. &
884	( t2m_daily(:) .GT. ncdgdd_temp(j)+ZeroCelsius ) )
885
886	gdd_midwinter(:,j) = &
887	gdd_midwinter(:,j) + &
888	dt * ( t2m_daily(:) - ( ncdgdd_temp(j)+ZeroCelsius ) )
889
890	ENDWHERE
891
892	ENDIF
893
894	ENDDO
895
896	!
897	!! 9. Calculate the number of chilling days (NCD) since leaves were lost (::ncd_dormance).
898	!! This variable is initialised to 0 at the beginning of the dormancy period (::time_lowgpp>0)
899	!! and increased by the stomate time step when the daily temperature is lower than the
900	!! PFT-dependent threshold ::ncdgdd_temp, which is defined for each PFT
901	!! in a table (::ncdgdd_temp_tab) in the module ::stomate_data. \n
902	!! It is set to "undef" when there is GPP (::time_lowgpp=0) as it is not needed during
903	!! the growing season. \n
904	!! The variable is used in the ::stomate_phenology module for the leaf onset model 'ncdgdd'.
905	!! Each PFT is looped over but the ::ncd_dormance is only calculated for those
906	!! PFTs for which a critical 'ncdgdd' temperature is defined, i.e. those which are
907	!! assigned to the 'ncdgdd' model.
908	!
909
910	DO j = 2,nvm ! Loop over # PFTs
911
912	IF ( ncdgdd_temp(j) .NE. undef ) THEN
913
914	!
915	!! 9.1 set to zero if undefined and there is no GPP
916	!
917
918	WHERE (gdd_init_date(:,1) .EQ. julian_diff)
919
920	ncd_dormance(:,j) = zero
921
922	ENDWHERE
923
924	!
925	!! 9.2 set to undef if there is GPP
926	!
927
928	WHERE ( when_growthinit(:,j) .EQ. zero )
929
930	ncd_dormance(:,j) = undef
931
932	ENDWHERE
933
934	!
935	!! 9.3 normal update, as described above, where ::ncd_dormance is defined
936	!
937
938	WHERE ( ( ncd_dormance(:,j) .NE. undef ) .AND. &
939	( t2m_daily(:) .LE. ncdgdd_temp(j)+ZeroCelsius ) )
940
941	ncd_dormance(:,j) = MIN( ncd_dormance(:,j) + dt, ncd_max )
942
943	ENDWHERE
944
945	ENDIF
946
947	ENDDO
948
949	!
950	!! 10. Calculate the number of growing days (NGD) since leaves were lost (::ngd_minus5).
951	!! This variable is initialised to 0 at the beginning of the dormancy period (::time_lowgpp>0)
952	!! and increased by the stomate time step when the daily temperature is higher than the threshold
953	!! -5 degrees C. \n
954	!! ::ngd_minus5 is further scaled by (::tau_ngd -dt / ::tau_ngd),
955	!! - ::tau_ngd is set to 50. in the module ::stomate_constants. \n
956	! Nicolas Viovy - not sure why this is...
957	!! The variable is used in the ::stomate_phenology module for the leaf onset model 'ngd'.
958	!! Each PFT is looped over.
959	!! If the NGD is less than zero, it is set to zero.
960	!
961
962	DO j = 2,nvm ! Loop over # PFTs
963
964	!
965	!! 10.1 Where there is GPP (i.e. ::time_lowgpp=0), set NGD to 0.
966	!! This means that we only take into account NGDs when the leaves are off
967	!
968
969	WHERE (gdd_init_date(:,1) .EQ. julian_diff)
970	ngd_minus5(:,j) = zero
971	ENDWHERE
972
973	!
974	!! 10.2 normal update, as described above.
975	!
976
977	WHERE ( t2m_daily(:) .GT. (ZeroCelsius - gdd_threshold) )
978	ngd_minus5(:,j) = ngd_minus5(:,j) + dt
979	ENDWHERE
980
981	ngd_minus5(:,j) = ngd_minus5(:,j) * ( tau_ngd - dt ) / tau_ngd
982
983	ENDDO
984
985	DO j = 2,nvm ! Loop over # PFTs
986	WHERE ( ABS(ngd_minus5(:,j)) .LT. EPSILON(zero) )
987	ngd_minus5(:,j) = zero
988	ENDWHERE
989	ENDDO
990
991	!
992	!! 11. Calculate the minimum humidity/relative soil moisture since dormance began (::hum_min_dormance)
993	!! and the time elapsed since this minimum (::time_hum_min). \n
994	!! The minimum moisture availability occurs when the monthly moisture availability, which is updated
995	!! daily earlier in this subroutine as previously described, is at a minimum during the dormancy period.
996	!! Therefore the ::hum_min_dormance is initialised to the monthly moisture availability
997	!! at the beginning of the dormancy period (i.e. if it was previously set to undefined and the
998	!! ::time_lowgpp>0) AND then whenever the monthly moisture availability is less than it has previously
999	!! been. \n
1000	!! Consequently, the time counter (::time_hum_min) is initialised to 0 at the beginning of the dormancy
1001	!! period (i.e. if it was previously set to undefined and the ::time_lowgpp>0) AND when the minimum
1002	!! moisture availability is reached, and is increased by the stomate time step every day throughout
1003	!! the dormancy period. \n
1004	!! ::time_hum_min is used in the ::stomate_phenology module for the leaf onset models 'moi' and 'moigdd'.
1005	!! Each PFT is looped over but the two variables are only calculated for those
1006	!! PFTs for which the critical parameter ::hum_min_time is defined, i.e.
1007	!! those which are assigned to the 'moi' or 'moigdd' models.
1008	!
1009
1010	DO j = 2,nvm ! Loop over # PFTs
1011
1012	IF ( hum_min_time(j) .NE. undef ) THEN
1013
1014	!
1015	!! 11.1 initialize if undefined and there is no GPP
1016	!
1017
1018	WHERE (when_growthinit(:,j) .EQ. zero)
1019
1020	time_hum_min(:,j) = zero
1021	hum_min_dormance(:,j) = moiavail_month(:,j)
1022
1023	ENDWHERE
1024
1025	!
1026	!! 11.3 normal update, as described above, where ::time_hum_min and ::hum_min_dormance are defined
1027	!
1028
1029	!! 11.3.1 increase time counter by the stomate time step
1030
1031	WHERE ( hum_min_dormance(:,j) .NE. undef )
1032
1033	time_hum_min(:,j) = time_hum_min(:,j) + dt
1034
1035	ENDWHERE
1036
1037	!! 11.3.2 set time counter to zero if minimum is reached
1038
1039	WHERE (( hum_min_dormance(:,j) .NE. undef ) .AND. &
1040	( moiavail_month(:,j) .LE. hum_min_dormance(:,j) ) )
1041
1042	hum_min_dormance(:,j) = moiavail_month(:,j)
1043	time_hum_min(:,j) = zero
1044
1045	ENDWHERE
1046
1047	ENDIF
1048
1049	ENDDO
1050
1051	!
1052	!! NOTE: Sections 12. to 14. compute slowly-varying, "long-term" (i.e. weekly/monthly)
1053	!! C fluxes (NPP, turnover, GPP) using the same linear relaxation method as in
1054	!! Sections 2. and 5. and described in the introduction to this section, as given
1055	!! by the following equation:
1056	!! \latexonly
1057	!! \input{season_lin_relax_eqn1.tex}
1058	!! \endlatexonly
1059	!! \n
1060	!! The following variables are calculated using the above equation, and the time constant
1061	!! is given for each.
1062	!
1063
1064	!
1065	!! 12. Update the "long term" NPP. (::npp_daily in gC/m^2/day, ::npp_longterm in gC/m^2/year.)
1066	!! The time constant is given by the parameter ::tau_longterm,
1067	!! which is set in ::stomate_data to be 3 * one year (in seconds, as
1068	!! described above). If the ::npp_longterm is less than zero then set it to zero. \n
1069	!! ::npp_longterm is used in ::stomate_lpj in the calculation of long-term
1070	!! vegetation dynamics and in ::stomate_lcchange in the calculation of
1071	!! land cover change. It is also used to calculate diagnose the hebivory activity in
1072	!! Section 22 of this subroutine.
1073	!
1074
1075	npp_longterm = ( npp_longterm * ( tau_longterm - dt ) + &
1076	(npp_dailyone_year) dt ) / &
1077	tau_longterm
1078
1079	DO j = 2,nvm ! Loop over # PFTs
1080	WHERE ( ABS(npp_longterm(:,j)) .LT. EPSILON(zero) )
1081	npp_longterm(:,j) = zero
1082	ENDWHERE
1083	ENDDO
1084
1085	!
1086	!! 13. Update the "long term" turnover rates (in gC/m^2/year).
1087	!! The time constant is given by the parameter ::tau_longterm,
1088	!! which is set in ::stomate_data to be 3 * one year (in seconds, as
1089	!! described above). If the ::turnover_longterm is less than zero then set it to zero.\n
1090	!! ::turnover_longterm is used in ::stomate_lpj and :: lpg_gap in the calculation
1091	!! of long-term vegetation dynamics.
1092	!
1093
1094	turnover_longterm(:,:,:,:) = ( turnover_longterm(:,:,:,:) * ( tau_longterm - dt ) + &
1095	(turnover_daily(:,:,:,:)one_year) dt ) / &
1096	tau_longterm
1097
1098	DO j = 2,nvm ! Loop over # PFTs
1099	WHERE ( ABS(turnover_longterm(:,j,:,:)) .LT. EPSILON(zero) )
1100	turnover_longterm(:,j,:,:) = zero
1101	ENDWHERE
1102	ENDDO
1103
1104	!
1105	!! 14. Update the "weekly" GPP (where there is vegetation), otherwise set to zero.
1106	!! The time constant is given by the parameter ::tau_gpp_week,
1107	!! which is set in ::stomate_data to be 7 days. If the ::gpp_week is
1108	!! less than zero then set it to zero. \n
1109	!! ::gpp_week is used to update the annual maximum weekly GPP (::maxgppweek_thisyear)
1110	!! in Section 16 of this subroutine, which is then used to update the variable
1111	!! ::maxgppweek_lastyear in Section 21 of this subroutine. Both ::gpp_week and
1112	!! ::maxgppweek_lastyear are used in Section 6 of this subroutine to calculate
1113	!! the onset and time-length of the dormancy period.
1114	! Note: Used to be weekly GPP divided by veget_max, i.e. per ground covered, but not anymore.
1115	!
1116
1117	WHERE ( veget_max .GT. zero )
1118
1119	gpp_week = ( gpp_week * ( tau_gpp_week - dt ) + &
1120	gpp_daily * dt ) / tau_gpp_week
1121
1122	ELSEWHERE
1123
1124	gpp_week = zero
1125
1126	ENDWHERE
1127
1128	DO j = 2,nvm ! Loop over # PFTs
1129	WHERE ( ABS(gpp_week(:,j)) .LT. EPSILON(zero) )
1130	gpp_week(:,j) = zero
1131	ENDWHERE
1132	ENDDO
1133
1134	!
1135	!! 15. Update the maximum and minimum moisture availabilities (::maxmoiavail_thisyear
1136	!! and ::minmoiavail_thisyear). If the daily moisture availability, ::moiavail_daily,
1137	!! which was calculated earlier in this subroutine, is greater or less than the current
1138	!! value of the maximum or minimum moisture availability, then set the value for the maximum
1139	!! or minimum to the current daily value respectively. \n
1140	!! ::maxmoiavail_thisyear and ::minmoiavail_thisyear are used to update the variables
1141	!! ::maxmoiavail_lastyear and ::minmoiavail_lastyear in Section 21 of this subroutine,
1142	!! which are used in the module ::stomate_phenology for the leaf onset models 'hum' and
1143	!! 'humgdd', and in the module ::stomate_turnover for the leaf senescence models 'dry'
1144	!! and 'mixed'.
1145	!
1146
1147	WHERE ( moiavail_daily .GT. maxmoiavail_thisyear )
1148	maxmoiavail_thisyear = moiavail_daily
1149	ENDWHERE
1150
1151	WHERE ( moiavail_daily .LT. minmoiavail_thisyear )
1152	minmoiavail_thisyear = moiavail_daily
1153	ENDWHERE
1154
1155	!
1156	!! 16. Update the annual maximum weekly GPP (::maxgppweek_thisyear), if the weekly GPP
1157	!! is greater than the current value of the annual maximum weekly GPP.
1158	!! The use of this variable is described in Section 14.
1159	!
1160
1161	WHERE ( gpp_week .GT. maxgppweek_thisyear )
1162	maxgppweek_thisyear = gpp_week
1163	ENDWHERE
1164
1165	!
1166	!! 17. Update the annual GDD0 by adding the current daily 2-meter temperature
1167	!! (multiplied by the stomate time step, which is one day), if it is greater
1168	!! than zero degrees C. \n
1169	!! This variable is mostly used in the module ::lpj_establish.
1170	!
1171
1172	WHERE ( t2m_daily .GT. ZeroCelsius )
1173	gdd0_thisyear = gdd0_thisyear + dt * ( t2m_daily - ZeroCelsius )
1174	ENDWHERE
1175
1176	!
1177	!! 18. Update the annual precipitation by adding the current daily precipitation
1178	!! amount (multiplied by the stomate time step, which is usually one day). \n
1179	!
1180
1181	precip_thisyear = precip_thisyear + dt * precip_daily
1182
1183
1184	!
1185	!! 19. Update the annual maximum leaf mass for each PFT (::lm_thisyearmax) and the maximum fractional
1186	!! plant cover if the LPJ DGVM is activated.
1187	!
1188
1189	!
1190	!! 19.1 If the LPJ DGVM is activated first the fraction of natural vegetation (::fracnat), i.e. non-
1191	!! agricultural vegetation (PFTs 2-11) is calculated for each PFT. Each PFT is looped over.
1192	!
1193	IF(ok_dgvm ) THEN
1194
1195	fracnat(:) = un
1196	DO j = 2,nvm ! Loop over # PFTs
1197	IF ( .NOT. natural(j) ) THEN
1198	fracnat(:) = fracnat(:) - veget_max(:,j)
1199	ENDIF
1200	ENDDO
1201
1202	ENDIF
1203
1204	!
1205	!! 19.2 If LPJ and STOMATE are activated, first the maximum fractional plant cover needs to be updated,
1206	!! and then this year's leaf biomass. Each PFT is looped over.
1207	!! Both are updated according to the linear relaxation method described above.
1208	!! \latexonly
1209	!! \input{season_lin_relax_eqn1.tex}
1210	!! \endlatexonly
1211	!! \n
1212	!! The time constant for this process is set as one year (in seconds) divided by the parameter
1213	!! ::leaflife_tab, which gives the leaf lifetime in years and is set for each PFT in the module
1214	!! ::stomate_constants. Each PFT is looped over. \n
1215	!
1216
1217	IF ( ok_stomate ) THEN
1218	IF(ok_dgvm ) THEN
1219	DO j=2,nvm ! Loop over # PFTs
1220
1221	!
1222	!! 19.2.1 Calculate maximum fractional plant cover (::maxfpc_lastyear).
1223	!! If natural vegetation is present in the grid cell, and the leaf
1224	!! biomass is greater than three-quarters of last year's leaf biomass, the maximum fractional plant
1225	!! cover for last year is updated. \n
1226	!! The short-term variable (Xs in the above equation) that is being used to update the long-term
1227	!! maximum fractional plant cover is the fractional cover of natural vegetation, specified as
1228	!! ::veget/::fracnat. Last year's value is then set to be this year's value.
1229	!
1230
1231	IF ( natural(j) .AND. ok_dgvm ) THEN
1232
1233	WHERE ( fracnat(:) .GT. min_stomate .AND. biomass(:,j,ileaf,icarbon).GT. lm_lastyearmax(:,j)*0.75 )
1234	maxfpc_lastyear(:,j) = ( maxfpc_lastyear(:,j) * ( one_year/leaflife_tab(j)- dt ) + &
1235	veget(:,j) / fracnat(:) * dt ) / (one_year/leaflife_tab(j))
1236	ENDWHERE
1237	maxfpc_thisyear(:,j) = maxfpc_lastyear(:,j) ! just to initialise value
1238
1239	ENDIF
1240
1241	!NV : correct initialization
1242	!!$ WHERE(biomass(:,j,ileaf,icarbon).GT. lm_lastyearmax(:,j)*0.75)
1243	!!$ lm_lastyearmax(:,j) = ( lm_lastyearmax(:,j) * ( one_year/leaflife_tab(j)- dt ) + &
1244	!!$ biomass(:,j,ileaf,icarbon) * dt ) / (one_year/leaflife_tab(j))
1245	!!$ ENDWHERE
1246	!!$ lm_thisyearmax(:,j)=lm_lastyearmax(:,j) ! just to initialise value
1247
1248
1249	!
1250	!! 19.2.2 Update this year's leaf biomass (::lm_thisyearmax).
1251	!! The short-term variable (Xs in the above equation) that is being used to update the long-term
1252	!! this year's leaf biomass is the leaf biomass pool (::biomass(i,j,ileaf).
1253	!
1254	WHERE (lm_thisyearmax(:,j) .GT. min_stomate)
1255	WHERE(biomass(:,j,ileaf,icarbon).GT. lm_thisyearmax(:,j)*0.75)
1256	lm_thisyearmax(:,j) = ( lm_thisyearmax(:,j) * ( one_year/leaflife_tab(j)- dt ) + &
1257	biomass(:,j,ileaf,icarbon) * dt ) / (one_year/leaflife_tab(j))
1258	ENDWHERE
1259	ELSEWHERE
1260	lm_thisyearmax(:,j) =biomass(:,j,ileaf,icarbon)
1261	ENDWHERE
1262
1263	ENDDO
1264
1265	ELSE
1266
1267	!
1268	!! 19.3 If LPJ DGVM is not activated but STOMATE is, the maximum leaf mass is set to be the same
1269	!! as the leaf biomass (::biomass(i,j,ileaf), without a change in the maximum fractional plant cover.
1270	!
1271	DO j = 2,nvm ! Loop over # PFTs
1272	WHERE ( biomass(:,j,ileaf,icarbon) .GT. lm_thisyearmax(:,j) )
1273	lm_thisyearmax(:,j) = biomass(:,j,ileaf,icarbon)
1274	ENDWHERE
1275	ENDDO
1276
1277	ENDIF !(ok_dgvm)
1278	ELSE
1279
1280	!
1281	!! 19.4 If STOMATE is not activated, the maximum leaf biomass is set to be the maximum possible
1282	!! LAI of the PFT.
1283	!
1284
1285	DO j = 2,nvm ! Loop over # PFTs
1286	!gmjc
1287	! lm_thisyearmax(:,j) = lai_max(j) / sla(j)
1288	lm_thisyearmax(:,j) = lai_max(j) / sla_calc(:,j)
1289	!end gmjc
1290	ENDDO
1291
1292	ENDIF !(ok_stomate)
1293
1294	!
1295	!! 20. Update the annual maximum fractional plant cover for each PFT if the current fractional cover
1296	!! (::veget), is larger than the current maximum value.
1297	!! ::veget is defined as fraction of total ground. Therefore, maxfpc_thisyear has the same unit.
1298	!
1299
1300	WHERE ( veget(:,:) .GT. maxfpc_thisyear(:,:) )
1301	maxfpc_thisyear(:,:) = veget(:,:)
1302	ENDWHERE
1303
1304	!
1305	!! 21. At the end of the every year, last year's maximum and minimum moisture availability,
1306	!! annual GDD0, annual precipitation, annual max weekly GPP, and maximum leaf mass are
1307	!! updated with the value calculated for the current year, and then their value reset.
1308	!! Either to zero for the maximum variables, or to ::large_value (set to be 1.E33 in ::
1309	!! the module ::stomate_constants) for the minimum variable ::minmoiavail_thisyear.
1310	!
1311
1312	IF ( EndOfYear ) THEN
1313
1314	!
1315	!! 21.1 Update last year's values. \n
1316	!! The variables ::maxmoiavail_lastyear, ::minmoiavail_lastyear and ::maxgppweek_lastyear
1317	!! are updated using the relaxation method:
1318	!! \latexonly
1319	!! \input{season_lin_relax_eqn1.tex}
1320	!! \endlatexonly
1321	!! \n
1322	!! where Xs is this year's value, dt (Delta-t) is set to 1 (year), and the time constant (tau) !??
1323	!! is set by the parameter ::tau_climatology, which is set to 20 years at the beginning
1324	!! of this subroutine. \n
1325	!! The other variables (::gdd0_lastyear, ::precip_lastyear, ::lm_lastyearmax and
1326	!! ::maxfpc_lastyear) are just replaced with this year's value.
1327	!
1328	!NVMODIF
1329	maxmoiavail_lastyear(:,:) = (maxmoiavail_lastyear(:,:)*(tau_climatology-1)+ maxmoiavail_thisyear(:,:))/tau_climatology
1330	minmoiavail_lastyear(:,:) = (minmoiavail_lastyear(:,:)*(tau_climatology-1)+ minmoiavail_thisyear(:,:))/tau_climatology
1331	maxgppweek_lastyear(:,:) =( maxgppweek_lastyear(:,:)*(tau_climatology-1)+ maxgppweek_thisyear(:,:))/tau_climatology
1332	! maxmoiavail_lastyear(:,:) = maxmoiavail_thisyear(:,:)
1333	! minmoiavail_lastyear(:,:) = minmoiavail_thisyear(:,:)
1334	! maxgppweek_lastyear(:,:) = maxgppweek_thisyear(:,:)
1335
1336	gdd0_lastyear(:) = gdd0_thisyear(:)
1337
1338	precip_lastyear(:) = precip_thisyear(:)
1339
1340	lm_lastyearmax(:,:) = lm_thisyearmax(:,:)
1341
1342	maxfpc_lastyear(:,:) = maxfpc_thisyear(:,:)
1343
1344	! Calculate Tseason
1345	Tseason(:) = zero
1346	WHERE ( Tseason_length(:) .GT. zero )
1347	Tseason(:) = Tseason_tmp(:) / Tseason_length(:)
1348	ENDWHERE
1349
1350	!
1351	!! 21.2 Reset new values for the "this year" variables. \n
1352	!! The maximum variables are set to zero and the minimum variable ::minmoiavail_thisyear
1353	!! is set to ::large_value (set to be 1.E33 in the module ::stomate_constants).
1354	!
1355
1356	maxmoiavail_thisyear(:,:) = zero
1357	minmoiavail_thisyear(:,:) = large_value
1358
1359	maxgppweek_thisyear(:,:) = zero
1360
1361	gdd0_thisyear(:) = zero
1362
1363	precip_thisyear(:) = zero
1364
1365	lm_thisyearmax(:,:) = zero
1366
1367	maxfpc_thisyear(:,:) = zero
1368
1369	Tseason_tmp(:) = zero
1370	Tseason_length(:) =zero
1371	Tmin_spring_time(:,:)=zero
1372	onset_date(:,:)=zero
1373
1374	!
1375	! 21.3 Special treatment for maxfpc. !??
1376	!! 21.3 Set the maximum fractional plant cover for non-natural vegetation
1377	!! (i.e. agricultural C3 and C4 - PFT 12 and 13) vegetation for last year to be zero.
1378	!
1379
1380	!
1381	! 21.3.1 Only take into account natural PFTs
1382	!
1383
1384	DO j = 2,nvm ! Loop over # PFTs
1385	IF ( .NOT. natural(j) ) THEN
1386	maxfpc_lastyear(:,j) = zero
1387	ENDIF
1388	ENDDO
1389
1390	! 21.3.2 In Stomate, veget is defined as a fraction of ground, not as a fraction
1391	! of total ground. maxfpc_lastyear will be compared to veget in lpj_light.
1392	! Therefore, we have to transform maxfpc_lastyear.
1393
1394
1395	! 21.3.3 The sum of the maxfpc_lastyear for natural PFT must not exceed fpc_crit (=.95).
1396	! However, it can slightly exceed this value as not all PFTs reach their maximum
1397	! fpc at the same time. Therefore, if sum(maxfpc_lastyear) for the natural PFTs
1398	! exceeds fpc_crit, we scale the values of maxfpc_lastyear so that the sum is
1399	! fpc_crit.
1400
1401	!!$ ! calculate the sum of maxfpc_lastyear
1402	!!$ sumfpc_nat(:) = zero
1403	!!$ DO j = 2,nvm ! Loop over # PFTs
1404	!!$ sumfpc_nat(:) = sumfpc_nat(:) + maxfpc_lastyear(:,j)
1405	!!$ ENDDO
1406	!!$
1407	!!$ ! scale so that the new sum is fpc_crit
1408	!!$ DO j = 2,nvm ! Loop over # PFTs
1409	!!$ WHERE ( sumfpc_nat(:) .GT. fpc_crit )
1410	!!$ maxfpc_lastyear(:,j) = maxfpc_lastyear(:,j) * (fpc_crit/sumfpc_nat(:))
1411	!!$ ENDWHERE
1412	!!$ ENDDO
1413
1414	ENDIF ! EndOfYear
1415
1416	!
1417	!! 22. Diagnose herbivore activity (::herbivores), determined as probability for a leaf to be
1418	!! eaten in a day (follows McNaughton et al., 1989). \n
1419	!! The amount of herbivore activity is used in the modules ::lpj_establish
1420	!! and ::stomate_turnover.
1421	!
1422
1423	!
1424	!! 22.1 First calculate the mean long-term leaf NPP in grid box, mean residence
1425	!! time (years) of green tissue to give the biomass available for herbivore consumption.
1426	!! Each PFT is looped over, though the calculation is only made for
1427	!! natural vegetation (PFTs 2-11).
1428	!
1429
1430	nlflong_nat(:,:) = zero
1431	weighttot(:,:) = zero
1432	green_age(:,:) = zero
1433	!
1434	DO j = 2,nvm ! Loop over # PFTs
1435	!
1436	IF ( natural(j) ) THEN
1437	!
1438	!! 22.1.1 Calculate the total weight of the leaves (::weighttot) as last year's leaf biomass
1439	!
1440	weighttot(:,j) = lm_lastyearmax(:,j)
1441
1442	!
1443	!! 22.1.2 Calculate the mean long-term leaf NPP as the long-term NPP calculated in Section 12 of
1444	!! this subroutine weighted by the leaf fraction (::leaf_frac), which is defined to be 0.33
1445	!! at the beginning of this subroutine.
1446	!
1447	nlflong_nat(:,j) = npp_longterm(:,j) * leaf_frac_hvc
1448
1449	!
1450	!! 22.1.3 Calculate the mean residence time of the green tissue (::green_age)
1451	!! This is calculated as the sum of 6 months
1452	!! for natural seasonal vegetation (i.e. PFTs 3, 6, 8-11), and 2 years for evergreen
1453	!! (PFTs 2, 4, 5, 8), multiplied by last year's leaf biomass for each PFT, divided by the
1454	!! total weight of the leaves for all PFTs.
1455	! This is a crude approximation. !!?? By whom?
1456	!! The difference between seasonal and evergreen vegetation is determined by the parameter
1457	!! ::pheno_model, which specifies the onset model of the PFT.
1458
1459	!
1460	IF ( pheno_model(j) .EQ. 'none' ) THEN
1461	green_age(:,j) = green_age_ever * lm_lastyearmax(:,j)
1462	ELSE
1463	green_age(:,j) = green_age_dec * lm_lastyearmax(:,j)
1464	ENDIF
1465	!
1466	ENDIF
1467	!
1468	ENDDO
1469	!
1470	WHERE ( weighttot(:,:) .GT. zero )
1471	green_age(:,:) = green_age(:,:) / weighttot(:,:)
1472	ELSEWHERE
1473	green_age(:,:) = un
1474	ENDWHERE
1475
1476	!
1477	!! 22.2 McNaughton et al. (1989) give herbivore consumption as a function of mean, long-term leaf NPP.
1478	!! as it gives an estimate of the edible biomass. The consumption of biomass by herbivores and the
1479	!! resultant herbivore activity are calculated following the equations:
1480	!! \latexonly
1481	!! \input{season_consumption_eqn4.tex}
1482	!! \endlatexonly
1483	!! \n
1484	!! and
1485	!! \latexonly
1486	!! \input{season_herbivore_eqn5.tex}
1487	!! \endlatexonly
1488	!! \n
1489	!
1490
1491	DO j = 2,nvm ! Loop over # PFTs
1492	!
1493	IF ( natural(j) ) THEN
1494	!
1495	WHERE ( nlflong_nat(:,j) .GT. zero )
1496	consumption(:) = hvc1 * nlflong_nat(:,j) ** hvc2
1497	herbivores(:,j) = one_year * green_age(:,j) * nlflong_nat(:,j) / consumption(:)
1498	ELSEWHERE
1499	herbivores(:,j) = 100000.
1500	ENDWHERE
1501	!
1502	ELSE
1503	!
1504	herbivores(:,j) = 100000.
1505	!
1506	ENDIF
1507	!
1508	ENDDO
1509	herbivores(:,ibare_sechiba) = zero
1510
1511	IF (printlev>=4) WRITE(numout,*) 'Leaving season'
1512
1513	END SUBROUTINE season
1514
1515	END MODULE stomate_season

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

Download in other formats: