Context Navigation

source: tags/ORCHIDEE_1_9_5_2/ORCHIDEE/src_stomate/stomate_season.f90 @ 8094

Last change on this file since 8094 was 248, checked in by martial.mancip, 13 years ago

Nicolas Viovy : Correction of 2 problems

the test on minimum annual water for grass establishment was removed
lm_max_lastyear was never reinitialised and then was not able to takes
into account for decrease of max lai for year to year

File size: 32.5 KB

Line
1	! Calculate long-term meteorological parameters from daily temperatures
2	! and precipitations (essentially for phenology)
3	!
4	! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/stomate_season.f90,v 1.15 2010/04/06 16:06:34 ssipsl Exp $
5	! IPSL (2006)
6	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
7	!
8	MODULE stomate_season
9
10	! modules used:
11
12	USE ioipsl
13	USE stomate_constants
14	USE constantes_veg
15
16	IMPLICIT NONE
17
18	! private & public routines
19
20	PRIVATE
21	PUBLIC season,season_clear
22
23	! first call
24	LOGICAL, SAVE :: firstcall = .TRUE.
25
26	CONTAINS
27
28
29	SUBROUTINE season_clear
30	firstcall=.TRUE.
31	END SUBROUTINE season_clear
32
33	SUBROUTINE season (npts, dt, EndOfYear, veget, veget_max, &
34	moiavail_daily, t2m_daily, tsoil_daily, soilhum_daily, &
35	precip_daily, npp_daily, biomass, turnover_daily, gpp_daily, when_growthinit, &
36	maxmoiavail_lastyear, maxmoiavail_thisyear, &
37	minmoiavail_lastyear, minmoiavail_thisyear, &
38	maxgppweek_lastyear, maxgppweek_thisyear, &
39	gdd0_lastyear, gdd0_thisyear, &
40	precip_lastyear, precip_thisyear, &
41	lm_lastyearmax, lm_thisyearmax, &
42	maxfpc_lastyear, maxfpc_thisyear, &
43	moiavail_month, moiavail_week, t2m_longterm, tlong_ref, t2m_month, t2m_week, &
44	tsoil_month, soilhum_month, &
45	npp_longterm, turnover_longterm, gpp_week, &
46	gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, time_lowgpp, &
47	time_hum_min, hum_min_dormance, herbivores)
48
49	!
50	! 0 declarations
51	!
52
53	! 0.1 input
54
55	! Domain size
56	INTEGER(i_std), INTENT(in) :: npts
57	! time step in days
58	REAL(r_std), INTENT(in) :: dt
59	! update yearly variables?
60	LOGICAL, INTENT(in) :: EndOfYear
61	! coverage fraction of a PFT. Here: fraction of total ground.
62	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget
63	! "maximal" coverage fraction of a PFT (for LAI -> infinity)
64	! Here: fraction of total ground.
65	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max
66	! Daily moisture availability
67	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: moiavail_daily
68	! Daily 2 meter temperature (K)
69	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_daily
70	! Daily soil temperature (K)
71	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: tsoil_daily
72	! Daily soil humidity
73	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: soilhum_daily
74	! Daily mean precipitation (mm/day)
75	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_daily
76	! daily net primary productivity (gC/m**2/day)
77	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: npp_daily
78	! biomass (gC/(m**2 of ground))
79	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: biomass
80	! Turnover rates (gC/m**2/day)
81	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: turnover_daily
82	! daily gross primary productivity (Here: gC/(m**2 of total ground)/day)
83	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: gpp_daily
84	! how many days ago was the beginning of the growing season
85	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: when_growthinit
86
87	! 0.2 modified fields
88
89	! last year's maximum moisture availability
90	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxmoiavail_lastyear
91	! this year's maximum moisture availability
92	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxmoiavail_thisyear
93	! last year's minimum moisture availability
94	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: minmoiavail_lastyear
95	! this year's minimum moisture availability
96	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: minmoiavail_thisyear
97	! last year's maximum weekly GPP
98	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxgppweek_lastyear
99	! this year's maximum weekly GPP
100	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxgppweek_thisyear
101	! last year's annual GDD0
102	REAL(r_std), DIMENSION(npts), INTENT(inout) :: gdd0_lastyear
103	! this year's annual GDD0
104	REAL(r_std), DIMENSION(npts), INTENT(inout) :: gdd0_thisyear
105	! last year's annual precipitation (mm/year)
106	REAL(r_std), DIMENSION(npts), INTENT(inout) :: precip_lastyear
107	! this year's annual precipitation (mm/year)
108	REAL(r_std), DIMENSION(npts), INTENT(inout) :: precip_thisyear
109	! last year's maximum leaf mass, for each PFT (gC/m**2)
110	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_lastyearmax
111	! this year's maximum leaf mass, for each PFT (gC/m**2)
112	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_thisyearmax
113	! last year's maximum fpc for each PFT, on ground
114	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxfpc_lastyear
115	! this year's maximum fpc for each PFT, on ground
116	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: maxfpc_thisyear
117	! "monthly" moisture availability
118	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_month
119	! "weekly" moisture availability
120	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_week
121	! "long term" 2-meter temperatures (K)
122	REAL(r_std), DIMENSION(npts), INTENT(inout) :: t2m_longterm
123	! "long term" refernce 2-meter temperatures (K)
124	REAL(r_std), DIMENSION(npts), INTENT(inout) :: tlong_ref
125	! "monthly" 2-meter temperatures (K)
126	REAL(r_std), DIMENSION(npts), INTENT(inout) :: t2m_month
127	! "weekly" 2-meter temperatures (K)
128	REAL(r_std), DIMENSION(npts), INTENT(inout) :: t2m_week
129	! "monthly" soil temperatures (K)
130	REAL(r_std), DIMENSION(npts,nbdl), INTENT(inout) :: tsoil_month
131	! "monthly" soil humidity
132	REAL(r_std), DIMENSION(npts,nbdl), INTENT(inout) :: soilhum_month
133	! "long term" net primary productivity (gC/m**2/year)
134	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: npp_longterm
135	! "long term" turnover rate (gC/m**2/year)
136	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: turnover_longterm
137	! "weekly" GPP (gC/day/(m**2 covered)
138	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gpp_week
139	! growing degree days, threshold -5 deg. C
140	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_m5_dormance
141	! growing degree days since midwinter
142	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_midwinter
143	! number of chilling days since leaves were lost
144	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ncd_dormance
145	! number of growing days
146	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ngd_minus5
147	! duration of dormance (d)
148	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: time_lowgpp
149	! time elapsed since strongest moisture availability (d)
150	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: time_hum_min
151	! minimum moisture during dormance
152	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: hum_min_dormance
153
154	! 0.3 output (diagnostic)
155
156	! time constant of probability of a leaf to be eaten by a herbivore (days)
157	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: herbivores
158
159	! 0.4 local
160
161	! indices
162	INTEGER(i_std) :: j
163	! rapport maximal GPP/GGP_max pour dormance
164	REAL(r_std), PARAMETER :: gppfrac_dormance = 0.2
165	!
166	!NVADD
167	! minimum gpp considered as not "lowgpp"
168	REAL(r_std), PARAMETER :: min_gpp_allowed = 0.3
169	! tau (year) for "climatologic variables
170	REAL(r_std), PARAMETER :: tau_climatology = 20
171	!ENDNVADD
172	! maximum ncd (d) (to avoid floating point underflows)
173	REAL(r_std) :: ncd_max
174	! parameters for herbivore activity
175	REAL(r_std), PARAMETER :: hvc1 = 0.019
176	REAL(r_std), PARAMETER :: hvc2 = 1.38
177	REAL(r_std), PARAMETER :: leaf_frac=.33
178	! sum of natural fpcs
179	REAL(r_std), DIMENSION(npts) :: sumfpc_nat
180	! weights
181	REAL(r_std), DIMENSION(npts) :: weighttot
182	! natural long-term leaf NPP ( gC/m**2/year)
183	REAL(r_std), DIMENSION(npts) :: nlflong_nat
184	! residence time of green tissue (years)
185	REAL(r_std), DIMENSION(npts) :: green_age
186	! herbivore consumption (gC/m**2/day)
187	REAL(r_std), DIMENSION(npts) :: consumption
188	! fraction of each gridcell occupied by natural vegetation
189	REAL(r_std), DIMENSION(npts) :: fracnat
190
191	! =========================================================================
192
193	IF (bavard.GE.3) WRITE(numout,*) 'Entering season'
194
195	!
196	! 1 Initializations
197	!
198	ncd_max = 3. * one_year
199
200	IF ( firstcall ) THEN
201
202	!
203	! 1.1 messages
204	!
205
206	IF ( bavard .GE. 1 ) THEN
207
208	WRITE(numout,*) 'season: '
209
210	WRITE(numout,*) ' > rapport maximal GPP/GGP_max pour dormance: ',gppfrac_dormance
211
212	WRITE(numout,*) ' > maximum possible ncd (d): ',ncd_max
213
214	WRITE(numout,*) ' > herbivore consumption C (gC/m2/day) as a function of NPP (gC/m2/d):'
215	WRITE(numout,) ' C=',hvc1,' NPP^',hvc2
216	WRITE(numout,) ' > for herbivores, suppose that ',leaf_frac100., &
217	'% of NPP is allocated to leaves'
218
219
220	ENDIF
221
222	!
223	! 1.2 Check whether longer-term meteorological parameters are initialized
224	! to zero
225	!
226
227	! 1.2.1 moisture availabilities
228
229	! 1.2.1.1 "monthly"
230	!MM PAS PARALLELISE!!
231	IF ( ABS( SUM( moiavail_month(:,2:nvm) ) ) .LT. min_stomate ) THEN
232
233	! in this case, set them it today's moisture availability
234	WRITE(numout,*) 'Warning! We have to initialize the ''monthly'' moisture availabilities.'
235	moiavail_month(:,:) = moiavail_daily(:,:)
236
237	ENDIF
238
239	! 1.2.1.2 "weekly"
240
241	IF ( ABS( SUM( moiavail_week(:,2:nvm) ) ) .LT. min_stomate ) THEN
242
243	! in this case, set them it today's moisture availability
244	WRITE(numout,*) 'Warning! We have to initialize the ''weekly'' moisture availabilities.'
245	moiavail_week(:,:) = moiavail_daily(:,:)
246
247	ENDIF
248
249	! 1.2.2 2-meter temperatures
250
251	! 1.2.2.1 "long term"
252
253	IF ( ABS( SUM( t2m_longterm(:) ) ) .LT. min_stomate ) THEN
254
255	! in this case, set them to today's temperature
256	WRITE(numout,*) 'Warning! We have to initialize the ''long term'' 2m temperatures.'
257	t2m_longterm(:) = t2m_daily(:)
258
259	ENDIF
260
261	! 1.2.2.2 "monthly"
262
263	IF ( ABS( SUM( t2m_month(:) ) ) .LT. min_stomate ) THEN
264
265	! in this case, set them to today's temperature
266	WRITE(numout,*) 'Warning! We have to initialize the ''monthly'' 2m temperatures.'
267	t2m_month(:) = t2m_daily(:)
268
269	ENDIF
270
271	! 1.2.2.3 "weekly"
272
273	IF ( ABS( SUM( t2m_week(:) ) ) .LT. min_stomate ) THEN
274
275	! in this case, set them to today's temperature
276	WRITE(numout,*) 'Warning! We have to initialize the ''weekly'' 2m temperatures.'
277	t2m_week(:) = t2m_daily(:)
278
279	ENDIF
280
281	! 1.2.3 "monthly" soil temperatures
282	!MM PAS PARALLELISE!!
283	IF ( ABS( SUM( tsoil_month(:,:) ) ) .LT. min_stomate ) THEN
284
285	! in this case, set them to today's temperature
286	WRITE(numout,*) 'Warning!'// &
287	' We have to initialize the ''monthly'' soil temperatures.'
288	tsoil_month(:,:) = tsoil_daily(:,:)
289
290	ENDIF
291
292	! 1.2.4 "monthly" soil humidity
293	IF ( ABS( SUM( soilhum_month(:,:) ) ) .LT. min_stomate ) THEN
294
295	! in this case, set them to today's humidity
296	WRITE(numout,*) 'Warning!'// &
297	' We have to initialize the ''monthly'' soil humidity.'
298	soilhum_month(:,:) = soilhum_daily(:,:)
299
300	ENDIF
301
302	! 1.2.5 growing degree days, threshold -5 deg C
303	IF ( ABS( SUM( gdd_m5_dormance(:,2:nvm) ) ) .LT. min_stomate ) THEN
304	WRITE(numout,*) 'Warning! Growing degree days (-5 deg) are initialized to ''undef''.'
305	gdd_m5_dormance(:,:) = undef
306	ENDIF
307
308	! 1.2.6 growing degree days since midwinter
309	IF ( ABS( SUM( gdd_midwinter(:,2:nvm) ) ) .LT. min_stomate ) THEN
310	WRITE(numout,*) 'Warning! Growing degree days since midwinter' // &
311	' are initialized to ''undef''.'
312	gdd_midwinter(:,:) = undef
313	ENDIF
314
315	! 1.2.7 number of chilling days since leaves were lost
316	IF ( ABS( SUM( ncd_dormance(:,2:nvm) ) ) .LT. min_stomate ) THEN
317	WRITE(numout,*) 'Warning! Number of chilling days is initialized to ''undef''.'
318	ncd_dormance(:,:) = undef
319	ENDIF
320
321	! 1.2.8 number of growing days, threshold -5 deg C
322	IF ( ABS( SUM( ngd_minus5(:,2:nvm) ) ) .LT. min_stomate ) THEN
323	WRITE(numout,*) 'Warning! Number of growing days (-5 deg) is initialized to 0.'
324	ENDIF
325
326	! 1.2.9 "long term" npp
327	IF ( ABS( SUM( npp_longterm(:,2:nvm) ) ) .LT. min_stomate ) THEN
328	WRITE(numout,*) 'Warning! Long term NPP is initialized to 0.'
329	ENDIF
330
331	! 1.2.10 "long term" turnover
332	IF ( ABS( SUM( turnover_longterm(:,2:nvm,:) ) ) .LT. min_stomate ) THEN
333	WRITE(numout,*) 'Warning! Long term turnover is initialized to 0.'
334	ENDIF
335
336	! 1.2.11 "weekly" GPP
337	IF ( ABS( SUM( gpp_week(:,2:nvm) ) ) .LT. min_stomate ) THEN
338	WRITE(numout,*) 'Warning! Weekly GPP is initialized to 0.'
339	ENDIF
340
341	! 1.2.12 minimum moisture availabilities
342	IF ( ABS( SUM( minmoiavail_thisyear(:,2:nvm) ) ) .LT. min_stomate ) THEN
343
344	! in this case, set them to a very high value
345	WRITE(numout,*) 'Warning! We have to initialize this year''s minimum '// &
346	'moisture availabilities.'
347	minmoiavail_thisyear(:,:) = large_value
348
349	ENDIF
350
351	!
352	! 1.3 reset flag
353	!
354
355	firstcall = .FALSE.
356
357	ENDIF
358
359	!
360	! 2 moisture availabilities
361	!
362
363	!
364	! 2.1 "monthly"
365	!
366
367	moiavail_month = ( moiavail_month * ( pheno_crit%tau_hum_month - dt ) + &
368	moiavail_daily * dt ) / pheno_crit%tau_hum_month
369
370	DO j = 2,nvm
371	WHERE ( ABS(moiavail_month(:,j)) .LT. EPSILON(zero) )
372	moiavail_month(:,j) = zero
373	ENDWHERE
374	ENDDO
375
376	!
377	! 2.2 "weekly"
378	!
379
380	moiavail_week = ( moiavail_week * ( pheno_crit%tau_hum_week - dt ) + &
381	moiavail_daily * dt ) / pheno_crit%tau_hum_week
382
383	DO j = 2,nvm
384	WHERE ( ABS(moiavail_week(:,j)) .LT. EPSILON(zero) )
385	moiavail_week(:,j) = zero
386	ENDWHERE
387	ENDDO
388
389	!
390	! 3 2-meter temperatures
391	!
392
393	!
394	! 3.1 "long term"
395	!
396
397	t2m_longterm = ( t2m_longterm * ( pheno_crit%tau_longterm - dt ) + &
398	t2m_daily * dt ) / pheno_crit%tau_longterm
399
400	WHERE ( ABS(t2m_longterm(:)) .LT. EPSILON(zero) )
401	t2m_longterm(:) = zero
402	ENDWHERE
403
404	CALL histwrite (hist_id_stomate, 'T2M_LONGTERM', itime, &
405	t2m_longterm, npts, hori_index)
406	!
407	! 3.2 "long term reference"
408	! This temperature is used for recalculating PFT-specific parameters such as
409	! critical photosynthesis temperatures of critical GDDs for phenology. This
410	! means that if the reference temperature varies, the PFTs adapt to them.
411	! Therefore the reference temperature can vary only if the vegetation is not
412	! static.
413	!
414
415	!!$ IF (control%ok_dgvm) THEN
416	tlong_ref(:) = MAX( tlong_ref_min, MIN( tlong_ref_max, t2m_longterm(:) ) )
417	!!$ ENDIF
418	!
419	! 3.3 "monthly"
420	!
421
422	t2m_month = ( t2m_month * ( pheno_crit%tau_t2m_month - dt ) + &
423	t2m_daily * dt ) / pheno_crit%tau_t2m_month
424
425	WHERE ( ABS(t2m_month(:)) .LT. EPSILON(zero) )
426	t2m_month(:) = zero
427	ENDWHERE
428
429	!
430	! 3.4 "weekly"
431	!
432
433	t2m_week = ( t2m_week * ( pheno_crit%tau_t2m_week - dt ) + &
434	t2m_daily * dt ) / pheno_crit%tau_t2m_week
435
436	WHERE ( ABS(t2m_week(:)) .LT. EPSILON(zero) )
437	t2m_week(:) = zero
438	ENDWHERE
439
440	!
441	! 4 ''monthly'' soil temperatures
442	!
443
444	tsoil_month = ( tsoil_month * ( pheno_crit%tau_tsoil_month - dt ) + &
445	tsoil_daily(:,:) * dt ) / pheno_crit%tau_tsoil_month
446
447	WHERE ( ABS(tsoil_month(:,:)) .LT. EPSILON(zero) )
448	tsoil_month(:,:) = zero
449	ENDWHERE
450
451	!
452	! 5 ''monthly'' soil humidity
453	!
454
455	soilhum_month = ( soilhum_month * ( pheno_crit%tau_soilhum_month - dt ) + &
456	soilhum_daily * dt ) / pheno_crit%tau_soilhum_month
457
458	WHERE ( ABS(soilhum_month(:,:)) .LT. EPSILON(zero) )
459	soilhum_month(:,:) = zero
460	ENDWHERE
461
462	!
463	! 6 dormance (d)
464	! when gpp is low, increase dormance time. Otherwise, set it to zero.
465	! NV: special case (3rd condition): plant is accumulating carbohydrates
466	! and does never use them. In this case, we allow the plant to
467	! detect a beginning of the growing season by declaring it dormant
468	!
469	!NVMODIF
470	DO j = 2,nvm
471	WHERE ( ( gpp_week(:,j) .LT. min_gpp_allowed ) .OR. &
472	( gpp_week(:,j) .LT. gppfrac_dormance * maxgppweek_lastyear(:,j) ) .OR. &
473	( ( when_growthinit(:,j) .GT. 2.*one_year ) .AND. &
474	( biomass(:,j,icarbres) .GT. biomass(:,j,ileaf)*4. ) ) )
475	! WHERE ( ( gpp_week(:,j) .EQ. zero ) .OR. &
476	! ( gpp_week(:,j) .LT. gppfrac_dormance * maxgppweek_lastyear(:,j) ) .OR. &
477	! ( ( when_growthinit(:,j) .GT. 2.*one_year ) .AND. &
478	! ( biomass(:,j,icarbres) .GT. biomass(:,j,ileaf)*4. ) ) )
479
480	time_lowgpp(:,j) = time_lowgpp(:,j) + dt
481
482	ELSEWHERE
483
484	time_lowgpp(:,j) = zero
485
486	ENDWHERE
487	ENDDO
488
489	!
490	! 7 growing degree days, threshold -5 deg C
491	!
492
493	DO j = 2,nvm
494
495	! only for PFTs for which critical gdd is defined
496	! gdd_m5_dormance is set to 0 at the end of the growing season. It is set to undef
497	! at the beginning of the growing season.
498
499	IF ( ALL(pheno_crit%gdd(j,:) .NE. undef) ) THEN
500
501	!
502	! 7.1 set to zero if undef and no gpp
503	!
504
505	WHERE ( ( time_lowgpp(:,j) .GT. zero ) .AND. ( gdd_m5_dormance(:,j) .EQ. undef ) )
506
507	gdd_m5_dormance(:,j) = zero
508
509	ENDWHERE
510
511	!
512	! 7.2 set to undef if there is gpp
513	!
514
515	WHERE ( time_lowgpp(:,j) .EQ. zero )
516
517	gdd_m5_dormance(:,j) = undef
518
519	ENDWHERE
520
521	!
522	! 7.3 normal update where gdd_m5_dormance is defined
523	!
524
525	WHERE ( ( t2m_daily(:) .GT. (ZeroCelsius-5.) ) .AND. &
526	( gdd_m5_dormance(:,j) .NE. undef ) )
527	gdd_m5_dormance(:,j) = gdd_m5_dormance(:,j) + &
528	dt * ( t2m_daily(:) - (ZeroCelsius-5.) )
529	ENDWHERE
530
531	WHERE ( gdd_m5_dormance(:,j) .NE. undef )
532	gdd_m5_dormance(:,j) = gdd_m5_dormance(:,j) * &
533	( pheno_crit%tau_gdd - dt ) / pheno_crit%tau_gdd
534	ENDWHERE
535
536	ENDIF
537
538	ENDDO
539
540	DO j = 2,nvm
541	WHERE ( ABS(gdd_m5_dormance(:,j)) .LT. EPSILON(zero) )
542	gdd_m5_dormance(:,j) = zero
543	ENDWHERE
544	ENDDO
545
546	!
547	! 8 growing degree days since midwinter
548	!
549
550	DO j = 2,nvm
551
552	! only for PFTs for which ncdgdd_crittemp is defined
553
554	IF ( pheno_crit%ncdgdd_temp(j) .NE. undef ) THEN
555
556	!
557	! 8.1 set to 0 if undef and if we detect "midwinter"
558	!
559
560	WHERE ( ( gdd_midwinter(:,j) .EQ. undef ) .AND. &
561	( t2m_month(:) .LT. t2m_week(:) ) .AND. &
562	( t2m_month(:) .LT. t2m_longterm(:) ) )
563
564	gdd_midwinter(:,j) = zero
565
566	ENDWHERE
567
568	!
569	! 8.2 set to undef if we detect "midsummer"
570	!
571
572	WHERE ( ( t2m_month(:) .GT. t2m_week(:) ) .AND. &
573	( t2m_month(:) .GT. t2m_longterm(:) ) )
574
575	gdd_midwinter(:,j) = undef
576
577	ENDWHERE
578
579	!
580	! 8.3 normal update
581	!
582
583	WHERE ( ( gdd_midwinter(:,j) .NE. undef ) .AND. &
584	( t2m_daily(:) .GT. pheno_crit%ncdgdd_temp(j)+ZeroCelsius ) )
585
586	gdd_midwinter(:,j) = &
587	gdd_midwinter(:,j) + &
588	dt * ( t2m_daily(:) - ( pheno_crit%ncdgdd_temp(j)+ZeroCelsius ) )
589
590	ENDWHERE
591
592	ENDIF
593
594	ENDDO
595
596	!
597	! 9 number of chilling days since leaves were lost
598	!
599
600	DO j = 2,nvm
601
602	IF ( pheno_crit%ncdgdd_temp(j) .NE. undef ) THEN
603
604	!
605	! 9.1 set to zero if undef and no gpp
606	!
607
608	WHERE ( ( time_lowgpp(:,j) .GT. zero ) .AND. ( ncd_dormance(:,j) .EQ. undef ) )
609
610	ncd_dormance(:,j) = zero
611
612	ENDWHERE
613
614	!
615	! 9.2 set to undef if there is gpp
616	!
617
618	WHERE ( time_lowgpp(:,j) .EQ. zero )
619
620	ncd_dormance(:,j) = undef
621
622	ENDWHERE
623
624	!
625	! 9.3 normal update where ncd_dormance is defined
626	!
627
628	WHERE ( ( ncd_dormance(:,j) .NE. undef ) .AND. &
629	( t2m_daily(:) .LE. pheno_crit%ncdgdd_temp(j)+ZeroCelsius ) )
630
631	ncd_dormance(:,j) = MIN( ncd_dormance(:,j) + dt, ncd_max )
632
633	ENDWHERE
634
635	ENDIF
636
637	ENDDO
638
639	!
640	! 10 number of growing days, threshold -5 deg C
641	!
642
643	DO j = 2,nvm
644
645	!
646	! 10.1 Where there is GPP, set ngd to 0
647	! This means that we only take into account ngds when the leaves are off
648	!
649
650	WHERE ( time_lowgpp(:,j) .LT. min_stomate )
651	ngd_minus5(:,j) = zero
652	ENDWHERE
653
654	!
655	! 10.2 normal update
656	!
657
658	WHERE ( t2m_daily(:) .GT. (ZeroCelsius-5.) )
659	ngd_minus5(:,j) = ngd_minus5(:,j) + dt
660	ENDWHERE
661
662	ngd_minus5(:,j) = ngd_minus5(:,j) * ( pheno_crit%tau_ngd - dt ) / pheno_crit%tau_ngd
663
664	ENDDO
665
666	DO j = 2,nvm
667	WHERE ( ABS(ngd_minus5(:,j)) .LT. EPSILON(zero) )
668	ngd_minus5(:,j) = zero
669	ENDWHERE
670	ENDDO
671	!
672	! 11 minimum humidity since dormance began and time elapsed since this minimum
673	!
674
675	DO j = 2,nvm
676
677	IF ( pheno_crit%hum_min_time(j) .NE. undef ) THEN
678
679	!
680	! 11.1 initialize if undef and no gpp
681	!
682
683	WHERE ( ( time_lowgpp(:,j) .GT. zero ) .AND. &
684	( ( time_hum_min(:,j) .EQ. undef ) .OR. ( hum_min_dormance(:,j) .EQ. undef ) ) )
685
686	time_hum_min(:,j) = zero
687	hum_min_dormance(:,j) = moiavail_month(:,j)
688
689	ENDWHERE
690
691	!
692	! 11.2 set to undef where there is gpp
693	!
694
695	WHERE ( time_lowgpp(:,j) .EQ. zero )
696
697	time_hum_min(:,j) = undef
698	hum_min_dormance(:,j) = undef
699
700	ENDWHERE
701
702	!
703	! 11.3 normal update where time_hum_min and hum_min_dormance are defined
704	!
705
706	! 11.3.1 increase time counter
707
708	WHERE ( ( time_hum_min(:,j) .NE. undef ) .AND. &
709	( hum_min_dormance(:,j) .NE. undef ) )
710
711	time_hum_min(:,j) = time_hum_min(:,j) + dt
712
713	ENDWHERE
714
715	! 11.3.2 set time to zero if minimum is reached
716
717	WHERE ( ( time_hum_min(:,j) .NE. undef ) .AND. &
718	( hum_min_dormance(:,j) .NE. undef ) .AND. &
719	( moiavail_month(:,j) .LE. hum_min_dormance(:,j) ) )
720
721	hum_min_dormance(:,j) = moiavail_month(:,j)
722	time_hum_min(:,j) = zero
723
724	ENDWHERE
725
726	ENDIF
727
728	ENDDO
729
730	!
731	! 12 "long term" NPP. npp_daily in gC/m2/day, npp_longterm in gC/m2/year.
732	!
733
734	npp_longterm = ( npp_longterm * ( pheno_crit%tau_longterm - dt ) + &
735	(npp_dailyone_year) dt ) / &
736	pheno_crit%tau_longterm
737
738	DO j = 2,nvm
739	WHERE ( ABS(npp_longterm(:,j)) .LT. EPSILON(zero) )
740	npp_longterm(:,j) = zero
741	ENDWHERE
742	ENDDO
743
744	!
745	! 13 "long term" turnover rates, in gC/m**2/year.
746	!
747
748	turnover_longterm = ( turnover_longterm * ( pheno_crit%tau_longterm - dt ) + &
749	(turnover_dailyone_year) dt ) / &
750	pheno_crit%tau_longterm
751
752	DO j = 2,nvm
753	WHERE ( ABS(turnover_longterm(:,j,:)) .LT. EPSILON(zero) )
754	turnover_longterm(:,j,:) = zero
755	ENDWHERE
756	ENDDO
757
758	!
759	! 14 "weekly" GPP, in gC/(m**2 covered)/day (!)
760	! i.e. divide daily gpp (in gC/m**2 of total ground/day) by vegetation fraction
761	! (m2 covered/m2 of total ground)
762	!
763
764	WHERE ( veget_max .GT. zero )
765
766	gpp_week = ( gpp_week * ( pheno_crit%tau_gpp_week - dt ) + &
767	gpp_daily * dt ) / pheno_crit%tau_gpp_week
768
769	ELSEWHERE
770
771	gpp_week = zero
772
773	ENDWHERE
774
775	DO j = 2,nvm
776	WHERE ( ABS(gpp_week(:,j)) .LT. EPSILON(zero) )
777	gpp_week(:,j) = zero
778	ENDWHERE
779	ENDDO
780
781	!
782	! 15 maximum and minimum moisture availabilities
783	!
784
785	WHERE ( moiavail_daily .GT. maxmoiavail_thisyear )
786	maxmoiavail_thisyear = moiavail_daily
787	ENDWHERE
788
789	WHERE ( moiavail_daily .LT. minmoiavail_thisyear )
790	minmoiavail_thisyear = moiavail_daily
791	ENDWHERE
792
793	!
794	! 16 annual maximum weekly GPP
795	!
796
797	WHERE ( gpp_week .GT. maxgppweek_thisyear )
798	maxgppweek_thisyear = gpp_week
799	ENDWHERE
800
801	!
802	! 17 annual GDD0
803	!
804
805	WHERE ( t2m_daily .GT. ZeroCelsius )
806	gdd0_thisyear = gdd0_thisyear + dt * ( t2m_daily - ZeroCelsius )
807	ENDWHERE
808
809	!
810	! 18 annual precipitation
811	!
812
813	precip_thisyear = precip_thisyear + dt * precip_daily
814
815	!
816	! 19 annual maximum leaf mass
817	! If STOMATE is not activated, this corresponds to the maximum possible
818	! LAI of the PFT
819	!
820
821	IF(control%ok_dgvm ) THEN
822
823	fracnat(:) = un
824	DO j = 2,nvm
825	IF ( .NOT. natural(j) ) THEN
826	fracnat(:) = fracnat(:) - veget_max(:,j)
827	ENDIF
828	ENDDO
829
830	ENDIF
831
832	IF ( control%ok_stomate ) THEN
833	IF(control%ok_dgvm ) THEN
834	DO j=2,nvm
835
836	IF ( natural(j) .AND. control%ok_dgvm ) THEN
837
838	WHERE ( fracnat(:) .GT. min_stomate .AND. biomass(:,j,ileaf).GT. lm_lastyearmax(:,j)*0.75 )
839	maxfpc_lastyear(:,j) = ( maxfpc_lastyear(:,j) * ( one_year/leaflife_tab(j)- dt ) + &
840	veget(:,j) / fracnat(:) * dt ) / (one_year/leaflife_tab(j))
841	ENDWHERE
842	maxfpc_thisyear(:,j) = maxfpc_lastyear(:,j) ! just to initialise value
843
844	ENDIF
845
846	!NV : correct initialization
847	!!$ WHERE(biomass(:,j,ileaf).GT. lm_lastyearmax(:,j)*0.75)
848	!!$ lm_lastyearmax(:,j) = ( lm_lastyearmax(:,j) * ( one_year/leaflife_tab(j)- dt ) + &
849	!!$ biomass(:,j,ileaf) * dt ) / (one_year/leaflife_tab(j))
850	!!$ ENDWHERE
851	!!$ lm_thisyearmax(:,j)=lm_lastyearmax(:,j) ! just to initialise value
852	WHERE (lm_thisyearmax(:,j) .GT. min_stomate)
853	WHERE(biomass(:,j,ileaf).GT. lm_thisyearmax(:,j)*0.75)
854	lm_thisyearmax(:,j) = ( lm_thisyearmax(:,j) * ( one_year/leaflife_tab(j)- dt ) + &
855	biomass(:,j,ileaf) * dt ) / (one_year/leaflife_tab(j))
856	ENDWHERE
857	ELSEWHERE
858	lm_thisyearmax(:,j) =biomass(:,j,ileaf)
859	ENDWHERE
860
861	ENDDO
862
863	ELSE
864
865	DO j = 2,nvm
866	WHERE ( biomass(:,j,ileaf) .GT. lm_thisyearmax(:,j) )
867	lm_thisyearmax(:,j) = biomass(:,j,ileaf)
868	ENDWHERE
869	ENDDO
870
871	ENDIF
872	ELSE
873
874	DO j = 2,nvm
875	lm_thisyearmax(:,j) = lai_max(j) / sla(j)
876	ENDDO
877
878	ENDIF
879
880	!
881	! 20 annual maximum fpc for each PFT
882	! "veget" is defined as fraction of total ground. Therefore, maxfpc_thisyear has
883	! the same unit.
884	!
885
886	WHERE ( veget(:,:) .GT. maxfpc_thisyear(:,:) )
887	maxfpc_thisyear(:,:) = veget(:,:)
888	ENDWHERE
889
890	!
891	! 21 Every year, replace last year's maximum and minimum moisture availability,
892	! annual GDD0, annual precipitation, annual max weekly GPP, and maximum leaf mass
893
894	IF ( EndOfYear ) THEN
895
896	!
897	! 21.1 replace old values
898	!
899	!NVMODIF
900	maxmoiavail_lastyear(:,:) = (maxmoiavail_lastyear(:,:)*(tau_climatology-1)+ maxmoiavail_thisyear(:,:))/tau_climatology
901	minmoiavail_lastyear(:,:) = (minmoiavail_lastyear(:,:)*(tau_climatology-1)+ minmoiavail_thisyear(:,:))/tau_climatology
902	maxgppweek_lastyear(:,:) =( maxgppweek_lastyear(:,:)*(tau_climatology-1)+ maxgppweek_thisyear(:,:))/tau_climatology
903	! maxmoiavail_lastyear(:,:) = maxmoiavail_thisyear(:,:)
904	! minmoiavail_lastyear(:,:) = minmoiavail_thisyear(:,:)
905	! maxgppweek_lastyear(:,:) = maxgppweek_thisyear(:,:)
906
907	gdd0_lastyear(:) = gdd0_thisyear(:)
908
909	precip_lastyear(:) = precip_thisyear(:)
910
911	lm_lastyearmax(:,:) = lm_thisyearmax(:,:)
912
913	maxfpc_lastyear(:,:) = maxfpc_thisyear(:,:)
914
915	!
916	! 21.2 reset new values
917	!
918
919	maxmoiavail_thisyear(:,:) = zero
920	minmoiavail_thisyear(:,:) = large_value
921
922	maxgppweek_thisyear(:,:) = zero
923
924	gdd0_thisyear(:) = zero
925
926	precip_thisyear(:) = zero
927
928	lm_thisyearmax(:,:) = zero
929
930	maxfpc_thisyear(:,:) = zero
931
932	!
933	! 21.3 Special treatment for maxfpc.
934	!
935
936	!
937	! 21.3.1 Only take into account natural PFTs
938	!
939
940	DO j = 2,nvm
941	IF ( .NOT. natural(j) ) THEN
942	maxfpc_lastyear(:,j) = zero
943	ENDIF
944	ENDDO
945
946	! 21.3.2 In Stomate, veget is defined as a fraction of ground, not as a fraction
947	! of total ground. maxfpc_lastyear will be compared to veget in lpj_light.
948	! Therefore, we have to transform maxfpc_lastyear.
949
950
951	! 21.3.3 The sum of the maxfpc_lastyear for natural PFT must not exceed fpc_crit (=.95).
952	! However, it can slightly exceed this value as not all PFTs reach their maximum
953	! fpc at the same time. Therefore, if sum(maxfpc_lastyear) for the natural PFTs
954	! exceeds fpc_crit, we scale the values of maxfpc_lastyear so that the sum is
955	! fpc_crit.
956
957	!!$ ! calculate the sum of maxfpc_lastyear
958	!!$ sumfpc_nat(:) = zero
959	!!$ DO j = 2,nvm
960	!!$ sumfpc_nat(:) = sumfpc_nat(:) + maxfpc_lastyear(:,j)
961	!!$ ENDDO
962	!!$
963	!!$ ! scale so that the new sum is fpc_crit
964	!!$ DO j = 2,nvm
965	!!$ WHERE ( sumfpc_nat(:) .GT. fpc_crit )
966	!!$ maxfpc_lastyear(:,j) = maxfpc_lastyear(:,j) * (fpc_crit/sumfpc_nat(:))
967	!!$ ENDWHERE
968	!!$ ENDDO
969
970	ENDIF ! EndOfYear
971
972	!
973	! 22 diagnose herbivore activity, determined through as probability for a leaf to be
974	! eaten in a day
975	! Follows McNaughton et al., Nat 341, 142-144, 1989.
976	!
977
978	!
979	! 22.1 first calculate mean long-term leaf NPP in grid box, mean residence
980	! time (years) of green tissue (i.e. tissue that will be eaten by
981	! herbivores) (crudely approximated: 6 months for seasonal and 2 years
982	! for evergreen) and mean length of growing season (6 months for
983	! seasonal and 1 year for evergreen).
984	!
985
986	nlflong_nat(:) = zero
987	weighttot(:) = zero
988	green_age(:) = zero
989	!
990	DO j = 2,nvm
991	!
992	IF ( natural(j) ) THEN
993	!
994	weighttot(:) = weighttot(:) + lm_lastyearmax(:,j)
995	nlflong_nat(:) = nlflong_nat(:) + npp_longterm(:,j) * leaf_frac
996	!
997	IF ( pheno_crit%pheno_model(j) .EQ. 'none' ) THEN
998	green_age(:) = green_age(:) + 2. * lm_lastyearmax(:,j)
999	ELSE
1000	green_age(:) = green_age(:) + .5 * lm_lastyearmax(:,j)
1001	ENDIF
1002	!
1003	ENDIF
1004	!
1005	ENDDO
1006	!
1007	WHERE ( weighttot(:) .GT. zero )
1008	green_age(:) = green_age(:) / weighttot(:)
1009	ELSEWHERE
1010	green_age(:) = un
1011	ENDWHERE
1012
1013	!
1014	! 22.2 McNaughton et al. give herbivore consumption as a function of annual leaf NPP.
1015	! The annual leaf NPP can give us an idea about the edible biomass:
1016	!
1017
1018	DO j = 2,nvm
1019	!
1020	IF ( natural(j) ) THEN
1021	!
1022	WHERE ( nlflong_nat(:) .GT. zero )
1023	consumption(:) = hvc1 * nlflong_nat(:) ** hvc2
1024	herbivores(:,j) = one_year * green_age(:) * nlflong_nat(:) / consumption(:)
1025	ELSEWHERE
1026	herbivores(:,j) = 100000.
1027	ENDWHERE
1028	!
1029	ELSE
1030	!
1031	herbivores(:,j) = 100000.
1032	!
1033	ENDIF
1034	!
1035	ENDDO
1036	herbivores(:,ibare_sechiba) = zero
1037
1038	IF (bavard.GE.4) WRITE(numout,*) 'Leaving season'
1039
1040	END SUBROUTINE season
1041
1042	END MODULE stomate_season

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

Download in other formats: