Context Navigation

source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_stomate/stomate_season.f90 @ 8398

Last change on this file since 8398 was 3642, checked in by josefine.ghattas, 8 years ago
Changes related to diagnostics according to the ORCHIDEE Retreat November 2016
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 67.9 KB

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

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

Download in other formats: