source: branches/publications/ORCHIDEE_CAN_r3069/src_stomate/lpj_pftinout.f90 @ 7475

Last change on this file since 7475 was 2442, checked in by sebastiaan.luyssaert, 10 years ago

DEV: replaced most STOPs by ipslerr. This version compiles but has not been tested

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 22.0 KB
Line 
1! =================================================================================================================================
2! MODULE       : lpj_pftinout
3!
4! CONTACT      : orchidee-help _at_ ipsl.jussieu.fr
5!
6! LICENCE      : IPSL (2006)
7! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8!
9!>\BRIEF       Introduce and eliminate PFT's from pixel
10!!
11!! \n DESCRIPTION: None
12!!
13!! RECENT CHANGE(S) : None
14!!
15!! REFERENCE(S) :
16!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
17!!        plant geography and terrestrial carbon cycling in the LPJ dynamic
18!!        global vegetation model, Global Change Biology, 9, 161-185.\n
19!! - Thonicke, K., S. Venevsky, et al. (2001), The role of fire disturbance
20!!        for global vegetation dynamics: coupling fire into a Dynamic Global
21!!        Vegetation Model, Global Ecology and Biogeography, 10, 661-677.\n
22!! - Haxeltine, A. and I. C. Prentice (1996), BIOME3: An equilibrium
23!!        terrestrial biosphere model based on ecophysiological constraints,
24!!        resource availability, and competition among plant functional types,
25!!        Global Biogeochemical Cycles, 10(4), 693-709.\n
26!! - Smith, B., I. C. Prentice, et al. (2001), Representation of vegetation
27!!        dynamics in the modelling of terrestrial ecosystems: comparing two
28!!        contrasting approaches within European climate space,
29!!        Global Ecology and Biogeography, 10, 621-637.\n
30!!
31!! SVN          :
32!! $HeadURL$
33!! $Date$
34!! $Revision$
35!! \n
36!_ ==============================================================================================================================
37
38MODULE lpj_pftinout
39
40  ! modules used:
41
42  USE ioipsl_para
43  USE stomate_data
44  USE pft_parameters
45  USE constantes
46
47  IMPLICIT NONE
48
49  ! private & public routines
50
51  PRIVATE
52  PUBLIC pftinout,pftinout_clear
53
54  LOGICAL, SAVE                       :: firstcall = .TRUE.                !! first call
55!$OMP THREADPRIVATE(firstcall)
56
57CONTAINS
58
59
60!! ================================================================================================================================
61!! SUBROUTINE  : pftinout_clear
62!!
63!>\BRIEF       Set flag ::firstcall to true and initialize the variables
64!_ ================================================================================================================================
65
66  SUBROUTINE pftinout_clear
67    firstcall = .TRUE.
68  END SUBROUTINE pftinout_clear
69
70!! ================================================================================================================================
71!! SUBROUTINE  : pftinout
72!!
73!>\BRIEF       Introduce and eliminate PFT's from pixel
74!!
75!! DESCRIPTION**3 : Introduction and elimination of PFTs on the basis of climate condition.
76!! For natural and woody PFTs the foliage projected coverage is calculated as follows:
77!! \latexonly
78!!  \input{equation_lpj_pftinout.tex}
79!! \endlatexonly
80!! \n
81!! where FPC is foliage projective cover (::fpc_nat), CN crown area (::cn_ind,
82!! @tex $ m^{2} $ @endtex), IND number of individuals (::ind, @tex $ m^{-2} $ @endtex,
83!! FRAC total fraction occupied by natural vegetation (::fracnat),
84!! @tex $ LM_{rm max} $ @endtex maximum leaf mass in last year
85!! (::lm_lastyearmax, @tex $ g C m^{-2} $ @endtex), SLA specific leaf area (sla,
86!! @tex $ m^{2} (g C)^{-1} $ @endtex), and coff coefficient (::ext_coeff). ::ext_coeff
87!! describes a property of the canopy (i.e. law of Lambert-Beer) and is defined in **2
88!!
89!! The foliage projective cover feeds into the calculation of the space available for
90!! expansion of existing and dispersion of new PFTs within a gridbox. In turn, available
91!! space is use to calculate the number of individuals with a PFT.
92!!
93!! Saplings are introduced under the condition that winter temperature is
94!! moderate, plant age is older than 1.25, (and for some PFTs at least one adjacent grid
95!! box exists for expansion), new saplings are introduced for narural PFT. In the simulation of
96!! agricultural grassland, if target PFT does not exist in the gridbox, it is introduced
97!! regardless of climate condition. When a new PFT is introduced CO_2 is taken from the
98!! atmosphere to account for CO_2 present in the seed and required by the germinated seeds
99!! to establish a sapling. These initial phases in ontology are not accounted for. However,
100!! by taking this small amount of CO2 from the atmosphere, mass balance closure for C is
101!! preserved.
102!!
103!! PFTs are eliminated under the condition that they are no longer adapted to the critical
104!! temperatures in winter. When a PFT is eliminated its number of indiviuals is set to zero and
105!! the rest of the elimination process is taken care of in lpj_kill.f90.
106!!
107!! RECENT CHANGE(S) : None
108!!
109!! MAIN OUTPUT VARIABLE(S): :: avail_tree (space availability for trees, unitless),   
110!! :: avail_grass (space availability for grasses, unitless), :: biomass (biomass, \f$gC m^{-2}\f$)
111!! and :: ind (density of individuals, \f$m^{-2}\f$)   
112!!
113!! REFERENCE(S) : None
114!!
115!! FLOWCHART    :
116!! \latexonly
117!!   \includegraphics[scale = 0.6]{pftinout.png}
118!! \endlatexonly
119!! \n
120!_ ================================================================================================================================
121
122
123
124
125
126
127  SUBROUTINE pftinout (npts, dt, adapted, regenerate, &
128       neighbours, veget_max, &
129       biomass, ind, cn_ind, age, leaf_frac, npp_longterm, lm_lastyearmax, senescence, &
130       PFTpresent, everywhere, when_growthinit, need_adjacent, RIP_time, &
131       co2_to_bm, &
132       avail_tree, avail_grass)
133   
134  !! 0. Variable and parameter declaration
135   
136    !! 0.1 Input variables
137   
138    INTEGER(i_std), INTENT(in)                                :: npts            !! Domain size - number of pixels (unitless)           
139    REAL(r_std), INTENT(in)                                   :: dt              !! Time step of vegetation dynamics for stomate
140                                                                                 !! (days)               
141    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: adapted         !! Winter not too cold (unitless)   
142    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: regenerate      !! Winter sufficiently cold (unitless) 
143    INTEGER(i_std), DIMENSION(npts,8), INTENT(in)             :: neighbours      !! Indices of the 8 neighbours of each grid point
144                                                                                 !! (unitless); 1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW,
145                                                                                 !! 7=W, 8=NW       
146    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: veget_max       !! "maximal" coverage fraction of a PFT (LAI ->
147                                                                                 !! infinity) on ground (unitless)       
148
149    !! 0.2 Output variables
150
151    REAL(r_std), DIMENSION(npts), INTENT(out)                 :: avail_tree      !! Space availability for trees (unitless)   
152    REAL(r_std), DIMENSION(npts), INTENT(out)                 :: avail_grass     !! Space availability for grasses (unitless)     
153
154
155    !! 0.3 Modified variables   
156
157    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass  !! Biomass (gC m^{-2})     
158    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: ind             !! Density of individuals (m^{-2})
159    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: cn_ind          !! Crown area of individuals (m^2)
160    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: age             !! Mean age (years)           
161    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac       !! Fraction of leaves in leaf age class (unitless) 
162    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: npp_longterm    !! "long term" net primary productivity
163                                                                                 !! (gC m^{-2} year^{-1})   
164    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: lm_lastyearmax  !! Last year's maximum leaf mass, for each PFT
165                                                                                 !! (gC m^{-2}) 
166    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)               :: senescence      !! Plant senescent for deciduous trees; .FALSE.
167                                                                                 !! if PFT is introduced or killed     
168    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)               :: PFTpresent      !! PFT exists   
169    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: everywhere      !! is the PFT everywhere in the grid box or very
170                                                                                 !! localized (unitless) 
171    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: when_growthinit !! how many days ago was the beginning of the
172                                                                                 !! growing season (days)
173    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)               :: need_adjacent   !! in order for this PFT to be introduced, does it
174                                                                                 !! have to be present in an adjacent grid box?   
175    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: RIP_time        !! How much time ago was the PFT eliminated for
176                                                                                 !! the last time (years)   
177    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: co2_to_bm       !! biomass uptaken (gC m^{-2} day^{-1})       
178 
179    !! 0.4 Local variables
180
181    REAL(r_std), DIMENSION(npts)                              :: avail           !! availability           
182    INTEGER(i_std)                                            :: i,j,m           !! indices     
183    REAL(r_std), DIMENSION(npts)                              :: sumfrac_wood    !! total woody vegetation cover     
184    INTEGER(i_std), DIMENSION(npts)                           :: n_present       !! number of adjacent grid cells where PFT is
185                                                                                 !! ubiquitous       
186    LOGICAL, DIMENSION(npts)                                  :: can_introduce   !! we can introduce this PFT   
187    REAL(r_std), DIMENSION(npts)                              :: fracnat         !! no real need for dimension(ntps) except for
188                                                                                 !! vectorisation         
189!_ ================================================================================================================================
190
191    IF (bavard.GE.3) WRITE(numout,*) 'Entering pftinout'
192
193  !! 1. Messages   
194
195    IF ( firstcall ) THEN
196
197       WRITE(numout,*) 'pftinout: Minimum space availability: ', min_avail
198
199       firstcall = .FALSE.
200
201    ENDIF
202
203  !! 2. Total woody fpc and space avaibility on grid
204
205    ! Only natural part of the grid cell\n
206    ! SZ bug correction MERGE: need to subtract agricultural area!
207    ! fraction of agricultural surface
208
209    !! 2.1 only natural PFT
210    fracnat(:) = un
211    DO j = 2,nvm ! Loop over # PFTs
212       IF ( .NOT. natural(j) ) THEN
213          fracnat(:) = fracnat(:) - veget_max(:,j)
214       ENDIF
215    ENDDO ! Loop over # PFTs
216
217    !! 2.2 Total woody fractional plant cover
218    sumfrac_wood(:) = zero
219    DO j = 2,nvm ! Loop over # PFTs
220       
221       ! SZ problem here: agriculture, not convinced that this representation of LPJ is correct
222       ! if agriculture is present, ind must be recalculated to correspond to the natural density...
223       ! since ind is per grid cell, can be achived by discounting for agricultura fraction
224       IF ( natural(j).AND.is_tree(j) ) THEN
225          WHERE(fracnat(:).GT.min_stomate) 
226                sumfrac_wood(:) = sumfrac_wood(:) + cn_ind(:,j) * ind(:,j) / fracnat(:) & 
227                     * ( un - exp( - lm_lastyearmax(:,j) * sla(j) * ext_coeff(j) ) )
228                !lai changed to lm_last
229          ENDWHERE
230       ENDIF
231    ENDDO ! Loop over # PFTs
232
233    !! 2.3 Space availability
234    avail_grass(:) = MAX( ( un - sumfrac_wood(:) ), min_avail )
235    avail_tree(:) = MAX( ( fpc_crit - sumfrac_wood(:) ), min_avail )
236
237  !! 3. Time since last elimination (y)
238
239    RIP_time = RIP_time + dt / one_year
240
241  !! 4. Agicultural PFTs
242
243    ! Agricultural PFTs are only present if they are prescribed
244    DO j = 2,nvm ! Loop over # PFTs
245       
246       IF ( .NOT. natural(j) ) THEN
247         
248          IF (bavard.GE.4) WRITE(numout,*) 'pftinout: Agricultural PFTs'
249         
250          !! 4.1 Agricultural trees
251          !      Agricultural trees are not treated for the moment
252          IF ( is_tree(j) ) THEN
253
254             CALL ipslerr_p (3,'lpj_pftinout',&
255                  'pftinout: Agricultural trees not treated.','','')
256
257          !! 4.2 Initialization of agricultural grass lands
258          !      Initialize parameter values of prescribed agricultural PFTs
259          ELSE
260
261             DO i = 1, npts ! Loop over # pixels - domain size
262
263                IF ( ( veget_max(i,j) .GT. min_stomate ) .AND. ( .NOT. PFTpresent(i,j) ) ) THEN
264
265                   ! prescribed, but not yet there.
266                   ind(i,j) = veget_max(i,j)
267                   biomass(i,j,:,:) = bm_sapl_old(j,:,:) * ind(i,j) /veget_max(i,j) 
268                   co2_to_bm(i,j) =  co2_to_bm(i,j) +SUM( biomass(i,j,:,icarbon) ) / dt
269                   PFTpresent(i,j) = .TRUE.
270                   everywhere(i,j) = un
271                   senescence(i,j) = .FALSE.
272                   age(i,j) = zero
273
274                ENDIF  ! prescribed, but PFT not yet present
275
276             ENDDO ! Loop over # pixels - domain size
277
278          ENDIF
279
280       ENDIF ! not natural
281
282    ENDDO ! Loop over # PFTs
283
284  !! 5 Eliminate PFTs
285
286    DO j = 2,nvm ! Loop over # PFTs
287
288       ! only for natural PFTs
289       IF ( natural(j) ) THEN
290
291          ! Number of individuals are set to zero in the condition of cold winter   
292          ! 'adapted_crit' critical value for being adapted = 1.-(1./euler); see 'stomate_constants.f90'
293          WHERE (  PFTpresent(:,j) .AND. ( adapted(:,j) .LT. adapted_crit ) )
294
295             ! PFT there, but not adapted any more (ex: winter too cold): kill
296             ! set number of individuals to zero - rest will be done in lpj_kill
297             ind(:,j) = zero
298
299          ENDWHERE
300                   WRITE(6,*) 'jfieozje lpj_pftinout',ind
301
302       ENDIF ! natural
303
304    ENDDO ! Loop over # PFTs
305
306  !! 6. Introduce PFTs
307
308    DO j = 2,nvm ! Loop over # PFTs
309
310       IF ( natural(j) ) THEN
311
312          ! space availability for this PFT
313          IF ( is_tree(j) ) THEN
314             avail(:) = avail_tree(:)
315          ELSE
316             avail(:) = avail_grass(:)
317          ENDIF
318         
319          !! 6.1 Check if PFT not present but (adapted and regenerative)     
320          can_introduce(:) = .FALSE.
321         
322          DO i = 1, npts ! Loop over # pixels - domain size
323
324             IF ( .NOT. PFTpresent(i,j) .AND. &
325                  ( adapted(i,j) .GT. adapted_crit ) .AND. &
326                  ( regenerate(i,j) .GT. regenerate_crit )  ) THEN
327
328                ! Seed are available nearby
329                IF ( need_adjacent(i,j) ) THEN
330
331                   !! 6.1.1 Climate allows introduction of the PFT but dispersion requires the
332                   !        presence of seeds. Seed are considered available if at least one
333                   !        neighbouring pixel is entirely invaded by the PFT. If that condition is
334                   !        satisfied, the PFT can establish in the new pixel.
335                   ! Count number of totally invaded neighbours.
336                   ! no loop so that it can vectorize
337                   n_present(i) = 0
338                   IF ( neighbours(i,1) .GT. 0 ) THEN
339                      IF ( everywhere(neighbours(i,1),j) .GE. un-min_stomate ) THEN
340                         n_present(i) = n_present(i)+1
341                      ENDIF
342                   ENDIF
343                   IF ( neighbours(i,3) .GT. 0 ) THEN
344                      IF ( everywhere(neighbours(i,3),j) .GE. un-min_stomate ) THEN
345                         n_present(i) = n_present(i)+1
346                      ENDIF
347                   ENDIF
348                   IF ( neighbours(i,5) .GT. 0 ) THEN
349                      IF ( everywhere(neighbours(i,5),j) .GE. un-min_stomate ) THEN
350                         n_present(i) = n_present(i)+1
351                      ENDIF
352                   ENDIF
353                   IF ( neighbours(i,7) .GT. 0 ) THEN
354                      IF ( everywhere(neighbours(i,7),j) .GE. un-min_stomate ) THEN
355                         n_present(i) = n_present(i)+1
356                      ENDIF
357                   ENDIF
358
359                   IF ( n_present(i) .GT. 0 ) THEN
360
361                      ! PFT is ubiquitous in at least one adjacent grid box
362                      can_introduce(i) = .TRUE.
363
364                   ENDIF
365
366                ELSE
367
368                   !! 6.1.2 No seed (trees) required for dispersion
369                   !        The PFT can establish without the presence of seed trees in
370                   !        neighbouring pixels.
371                   can_introduce(i) = .TRUE.
372
373                ENDIF ! do we have to look at the neighbours?
374
375             ENDIF ! we'd like to introduce the PFT
376
377          ENDDO ! Loop over # pixels - domain size
378
379          !! 6.2 Has the PFT been eliminated lately?
380          !      Additional test whether the PFT has been eliminated lately, i.e.
381          !      less than 1.25 years ago. Do not only take full years as success of
382          !      introduction, as introduction might depend on season.
383          WHERE ( RIP_time(:,j) .LT. RIP_time_min )
384
385             ! PFT was eliminated lately - cannot reintroduce
386             can_introduce(:) = .FALSE.
387
388          ENDWHERE
389
390          !! 6.3 Introduce that PFT where possible
391          !      "can_introduce" means that it either exists in neighbouring grid boxes
392          !      or that we do not look at neighbours, that it has not been eliminated
393          !      lately, and, of course, that the climate is good for that PFT.
394          WHERE ( can_introduce(:) )
395             
396             PFTpresent(:,j) = .TRUE.
397             
398             senescence(:,j) = .FALSE.
399             
400             ! introduce at least a few saplings, even if canopy is closed
401             ! initial density of individuals (ind_0) = 0.02, see 'stomate_constant.f90'
402             ind(:,j) = ind_0 * (dt/one_year) * avail(:)
403             
404             WHERE(veget_max(:,j) .GT. min_stomate)
405                biomass(:,j,ileaf,icarbon) = bm_sapl_old(j,ileaf,icarbon) * ind(:,j) /veget_max(:,j)
406                biomass(:,j,isapabove,icarbon) = bm_sapl_old(j,isapabove,icarbon) * ind(:,j) /veget_max(:,j)
407                biomass(:,j,isapbelow,icarbon) = bm_sapl_old(j,isapbelow,icarbon) * ind(:,j)/veget_max(:,j)
408                biomass(:,j,iheartabove,icarbon) = bm_sapl_old(j,iheartabove,icarbon) * ind(:,j)/veget_max(:,j)
409                biomass(:,j,iheartbelow,icarbon) = bm_sapl_old(j,iheartbelow,icarbon) * ind(:,j)/veget_max(:,j)
410                biomass(:,j,iroot,icarbon) = bm_sapl_old(j,iroot,icarbon) * ind(:,j)/veget_max(:,j)
411                biomass(:,j,ifruit,icarbon) = bm_sapl_old(j,ifruit,icarbon) * ind(:,j)/veget_max(:,j)
412                biomass(:,j,icarbres,icarbon) = bm_sapl_old(j,icarbres,icarbon) * ind(:,j)/veget_max(:,j)
413             ELSEWHERE             
414                biomass(:,j,ileaf,icarbon) = bm_sapl_old(j,ileaf,icarbon) * ind(:,j)
415                biomass(:,j,isapabove,icarbon) = bm_sapl_old(j,isapabove,icarbon) * ind(:,j)
416                biomass(:,j,isapbelow,icarbon) = bm_sapl_old(j,isapbelow,icarbon) * ind(:,j)
417                biomass(:,j,iheartabove,icarbon) = bm_sapl_old(j,iheartabove,icarbon) * ind(:,j)
418                biomass(:,j,iheartbelow,icarbon) = bm_sapl_old(j,iheartbelow,icarbon) * ind(:,j)
419                biomass(:,j,iroot,icarbon) = bm_sapl_old(j,iroot,icarbon) * ind(:,j)
420                biomass(:,j,ifruit,icarbon) = bm_sapl_old(j,ifruit,icarbon) * ind(:,j)
421                biomass(:,j,icarbres,icarbon) = bm_sapl_old(j,icarbres,icarbon) * ind(:,j)
422             END WHERE
423           
424             co2_to_bm(:,j) = &
425                  co2_to_bm(:,j) / dt * &
426                  ( biomass(:,j,ileaf,icarbon) + biomass(:,j,isapabove,icarbon) + &
427                  biomass(:,j,isapbelow,icarbon) + biomass(:,j,iheartabove,icarbon) + &
428                  biomass(:,j,iheartbelow,icarbon) + biomass(:,j,iroot,icarbon) + &
429                  biomass(:,j,ifruit,icarbon) + biomass(:,j,icarbres,icarbon) )
430
431             when_growthinit(:,j) = large_value
432
433             age(:,j) = zero
434
435             ! all leaves are young
436             leaf_frac(:,j,1) = un
437
438             ! non-zero "long term" npp and last year's leaf mass for saplings -
439             ! so they won't be killed off by gap or kill
440             npp_longterm(:,j) = npp_longterm_init
441
442             lm_lastyearmax(:,j) = bm_sapl_old(j,ileaf,icarbon) * ind(:,j)
443
444          ENDWHERE    ! we can introduce the PFT
445
446          !! 6.4 Expansion of the PFT within the grid box
447          !      PFT expansion/dispersion to a new grid box should not be confused with
448          !      expansion in areal coverage
449          IF ( treat_expansion ) THEN
450
451             WHERE ( can_introduce(:) )
452
453                ! low value at the beginning
454                everywhere(:,j) = everywhere_init
455             ENDWHERE
456
457          ELSE
458
459             ! expansion is not treated
460             WHERE ( can_introduce(:) )
461                everywhere(:,j) = un
462             ENDWHERE
463
464          ENDIF ! treat expansion
465
466       ENDIF ! only natural PFTs
467
468    ENDDO ! Loop over # PFTs
469
470  !! 7. If a PFT has been present once in a grid box, we suppose that it will survive
471
472    !   If a PFT has been present once in a grid box, we suppose that it will survive
473    !   in isolated places (e.g., an oasis) within that grid box, even if it gets
474    !   officially eliminated from it later. That means that if climate becomes favorable
475    !   again, it will not need to get seeds from adjacent grid cells.
476    WHERE ( PFTpresent )
477       need_adjacent = .FALSE.
478    ENDWHERE
479
480    IF (bavard.GE.4) WRITE(numout,*) 'Leaving pftinout'
481
482  END SUBROUTINE pftinout
483
484END MODULE lpj_pftinout
Note: See TracBrowser for help on using the repository browser.