Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-261/src_stomate/stomate_season.f90 @ 8692

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

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

Download in other formats: