Context Navigation

stomate_season.f90 @ 7346

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

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

Context Navigation

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

Download in other formats: