Context Navigation

stomate_season.f90 @ 8367

Last change on this file since 8367 was 7215, checked in by nicolas.vuichard, 3 years ago
correction of the last commit. No change in stomate_season to consider
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 73.6 KB

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

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

Context Navigation

source: branches/ORCHIDEE_3_CMIP6/ORCHIDEE/src_stomate/stomate_season.f90 @ 8367

Download in other formats: