! establishment routine ! Suppose seed pool >> establishment rate. ! ! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/lpj_establish.f90,v 1.9 2009/01/06 15:01:25 ssipsl Exp $ ! IPSL (2006) ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC ! MODULE lpj_establish ! modules used: USE ioipsl USE stomate_constants USE constantes_veg IMPLICIT NONE ! private & public routines PRIVATE PUBLIC establish,establish_clear ! first call LOGICAL, SAVE :: firstcall = .TRUE. CONTAINS SUBROUTINE establish_clear firstcall = .TRUE. END SUBROUTINE establish_clear SUBROUTINE establish (npts, dt, PFTpresent, regenerate, & neighbours, resolution, need_adjacent, herbivores, & precip_annual, gdd0, lm_lastyearmax, & cn_ind, lai, avail_tree, avail_grass, & leaf_age, leaf_frac, & ind, biomass, age, everywhere, co2_to_bm,veget_max) ! ! 0 declarations ! ! 0.1 input ! Domain size INTEGER(i_std), INTENT(in) :: npts ! Time step of vegetation dynamics (days) REAL(r_std), INTENT(in) :: dt ! Is pft there LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: PFTpresent ! Winter sufficiently cold REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: regenerate ! indices of the 8 neighbours of each grid point (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) INTEGER(i_std), DIMENSION(npts,8), INTENT(in) :: neighbours ! resolution at each grid point in m (1=E-W, 2=N-S) REAL(r_std), DIMENSION(npts,2), INTENT(in) :: resolution ! in order for this PFT to be introduced, does it have to be present in an ! adjacent grid box? LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: need_adjacent ! time constant of probability of a leaf to be eaten by a herbivore (days) REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: herbivores ! annual precipitation (mm/year) REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_annual ! growing degree days (C) REAL(r_std), DIMENSION(npts), INTENT(in) :: gdd0 ! last year's maximum leaf mass, for each PFT (gC/(m**2 of ground)) REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lm_lastyearmax ! crown area of individuals (m**2) REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: cn_ind ! leaf area index OF AN INDIVIDUAL PLANT REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lai ! space availability for trees REAL(r_std), DIMENSION(npts), INTENT(in) :: avail_tree ! space availability for grasses REAL(r_std), DIMENSION(npts), INTENT(in) :: avail_grass ! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max ! 0.2 modified fields ! leaf age (days) REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_age ! fraction of leaves in leaf age class REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac ! Number of individuals / m2 REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind ! biomass (gC/(m**2 of ground)) REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: biomass ! mean age (years) REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age ! is the PFT everywhere in the grid box or very localized (after its introduction) REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere ! biomass uptaken (gC/(m**2 of total ground)/day) !NV passage 2D REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_to_bm ! 0.3 local ! time during which a sapling can be entirely eaten by herbivores (d) REAL(r_std) :: tau_eatup ! new fpc ( foliage protected cover: fractional coverage ) REAL(r_std), DIMENSION(npts,nvm) :: fpc_nat ! maximum tree establishment rate, based on climate only REAL(r_std), DIMENSION(npts) :: estab_rate_max_climate_tree ! maximum grass establishment rate, based on climate only REAL(r_std), DIMENSION(npts) :: estab_rate_max_climate_grass ! maximum tree establishment rate, based on climate and fpc REAL(r_std), DIMENSION(npts) :: estab_rate_max_tree ! maximum grass establishment rate, based on climate and fpc REAL(r_std), DIMENSION(npts) :: estab_rate_max_grass ! total natural fpc REAL(r_std), DIMENSION(npts) :: sumfpc ! total woody fpc REAL(r_std), DIMENSION(npts) :: sumfpc_wood ! for trees, measures the total concurrence for available space REAL(r_std), DIMENSION(npts) :: spacefight_tree ! for grasses, measures the total concurrence for available space REAL(r_std), DIMENSION(npts) :: spacefight_grass ! change in number of individuals /m2 per time step (per day in history file) REAL(r_std), DIMENSION(npts,nvm) :: d_ind ! biomass increase (gC/(m**2 of ground)) REAL(r_std), DIMENSION(npts) :: bm_new ! stem diameter (m) REAL(r_std), DIMENSION(npts) :: dia ! temporary variable REAL(r_std), DIMENSION(npts) :: b1 ! new sap mass (gC/(m**2 of ground)) REAL(r_std), DIMENSION(npts) :: sm2 ! woodmass of an individual REAL(r_std), DIMENSION(npts) :: woodmass ! ratio of hw(above) to total hw, sm(above) to total sm REAL(r_std), DIMENSION(npts) :: sm_at ! reduction factor for establishment if many trees or grasses are present REAL(r_std), DIMENSION(npts) :: factor ! from how many sides is the grid box invaded INTEGER(i_std) :: nfrontx INTEGER(i_std) :: nfronty ! daily establishment rate is large compared to present number of individuals LOGICAL, DIMENSION(npts) :: many_new ! indices INTEGER(i_std) :: i,j,k,m ! ========================================================================= IF (bavard.GE.3) WRITE(numout,*) 'Entering establish' ! ! 1 messages and initialization ! tau_eatup = one_year/2. IF ( firstcall ) THEN WRITE(numout,*) 'establish:' WRITE(numout,*) ' > time during which a sapling can be entirely eaten by herbivores (d): ', & tau_eatup firstcall = .FALSE. ENDIF ! ! 2 recalculate fpc ! ! ! 2.1 Only natural part of the grid cell ! DO j = 2,nvm IF ( natural(j) ) THEN DO i = 1, npts IF (lai(i,j) == val_exp) THEN fpc_nat(i,j) = cn_ind(i,j) * ind(i,j) ELSE fpc_nat(i,j) = cn_ind(i,j) * ind(i,j) * ( un - exp( -lai(i,j) * ext_coeff(j) ) ) ENDIF ENDDO ELSE fpc_nat(:,j) = zero ENDIF ENDDO ! ! 2.2 total natural fpc on grid ! sumfpc(:) = SUM( fpc_nat(:,:), DIM=2 ) ! ! 2.3 total woody fpc on grid and number of regenerative tree pfts ! sumfpc_wood(:) = zero spacefight_tree(:) = zero DO j = 2,nvm IF ( tree(j) .AND. natural(j) ) THEN ! total woody fpc WHERE ( PFTpresent(:,j) ) sumfpc_wood(:) = sumfpc_wood(:) + fpc_nat(:,j) ENDWHERE ! how many trees are competing? Count a PFT fully only if it is present ! on the whole grid box. WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) ) spacefight_tree(:) = spacefight_tree(:) + everywhere(:,j) ENDWHERE ENDIF ENDDO ! ! 2.4 number of natural grasses ! spacefight_grass(:) = zero DO j = 2,nvm IF ( .NOT. tree(j) .AND. natural(j) ) THEN ! how many grasses are competing? Count a PFT fully only if it is present ! on the whole grid box. WHERE ( PFTpresent(:,j) ) spacefight_grass(:) = spacefight_grass(:) + everywhere(:,j) ENDWHERE ENDIF ENDDO ! ! 3 establishment rate ! ! ! 3.1 maximum establishment rate, based on climate only ! WHERE ( ( precip_annual(:) .GE. precip_crit ) .AND. ( gdd0(:) .GE. gdd_crit ) ) estab_rate_max_climate_tree(:) = estab_max_tree estab_rate_max_climate_grass(:) = estab_max_grass ELSEWHERE estab_rate_max_climate_tree(:) = zero estab_rate_max_climate_grass(:) = zero ENDWHERE ! ! 3.2 reduce maximum tree establishment rate if many trees present. ! In the original DGVM, this is done using a step function which yields a ! reduction by factor 4 if sumfpc_wood(i) .GT. fpc_crit - 0.05. ! This can lead to small oscillations (without consequences however). ! Here, a steady linear transition is used between fpc_crit-0.075 and ! fpc_crit-0.025. ! factor(:) = un - 15. * ( sumfpc_wood(:) - (fpc_crit-.075) ) factor(:) = MAX( 0.25_r_std, MIN( 1._r_std, factor(:) ) ) estab_rate_max_tree(:) = estab_rate_max_climate_tree(:) * factor(:) ! ! 3.3 Modulate grass establishment rate. ! If canopy is not closed (fpc < fpc_crit-0.05), normal establishment. ! If canopy is closed, establishment is reduced by a factor 4. ! Factor is linear between these two bounds. ! This is different from the original DGVM where a step function is ! used at fpc_crit-0.05 (This can lead to small oscillations, ! without consequences however). ! factor(:) = un - 15. * ( sumfpc(:) - (fpc_crit-.05) ) factor(:) = MAX( 0.25_r_std, MIN( 1._r_std, factor(:) ) ) estab_rate_max_grass(:) = estab_rate_max_climate_grass(:) * factor(:) ! ! 4 do establishment for natural PFTs ! d_ind(:,:) = zero DO j = 2,nvm ! only for natural PFTs IF ( natural(j) ) THEN ! ! 4.1 PFT expansion across the grid box. Not to be confused with areal ! coverage. ! IF ( treat_expansion ) THEN ! only treat plants that are regenerative and present and still can expand DO i = 1, npts IF ( PFTpresent(i,j) .AND. & ( everywhere(i,j) .LT. un ) .AND. & ( regenerate(i,j) .GT. regenerate_crit ) ) THEN ! from how many sides is the grid box invaded (separate x and y directions ! because resolution may be strongly anisotropic) ! ! For the moment we only look into 4 direction but that can be extanded (JP) ! nfrontx = 0 IF ( neighbours(i,3) .GT. 0 ) THEN IF ( everywhere(neighbours(i,3),j) .GT. 1.-min_stomate ) nfrontx = nfrontx+1 ENDIF IF ( neighbours(i,7) .GT. 0 ) THEN IF ( everywhere(neighbours(i,7),j) .GT. 1.-min_stomate ) nfrontx = nfrontx+1 ENDIF nfronty = 0 IF ( neighbours(i,1) .GT. 0 ) THEN IF ( everywhere(neighbours(i,1),j) .GT. 1.-min_stomate ) nfronty = nfronty+1 ENDIF IF ( neighbours(i,5) .GT. 0 ) THEN IF ( everywhere(neighbours(i,5),j) .GT. 1.-min_stomate ) nfronty = nfronty+1 ENDIF everywhere(i,j) = & everywhere(i,j) + migrate(j) * dt/one_year * & ( nfrontx / resolution(i,1) + nfronty / resolution(i,2) ) IF ( .NOT. need_adjacent(i,j) ) THEN ! in that case, we also assume that the PFT expands from places within ! the grid box (e.g., oasis). everywhere(i,j) = & everywhere(i,j) + migrate(j) * dt/one_year * & 2. * SQRT( pi*everywhere(i,j)/(resolution(i,1)*resolution(i,2)) ) ENDIF everywhere(i,j) = MIN( everywhere(i,j), 1._r_std ) ENDIF ENDDO ENDIF ! treat expansion? ! ! 4.2 establishment rate ! - Is lower if the PFT is only present in a small part of the grid box ! (after its introduction), therefore multiplied by "everywhere". ! - Is divided by the number of PFTs that compete ("spacefight"). ! - Is modulated by space availability (avail_tree, avail_grass). ! IF ( tree(j) ) THEN ! 4.2.1 present and regenerative trees WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) ) d_ind(:,j) = estab_rate_max_tree(:)*everywhere(:,j)/spacefight_tree(:) * & avail_tree(:) * dt/one_year ENDWHERE ELSE ! 4.2.2 present and regenerative grasses WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) ) d_ind(:,j) = estab_rate_max_grass(:)*everywhere(:,j)/spacefight_grass(:) * & avail_grass(:) * dt/one_year ENDWHERE ENDIF ! tree/grass ! ! 4.3 herbivores reduce establishment rate ! We suppose that saplings are vulnerable during a given time after establishment. ! This is taken into account by preventively reducing the establishment rate. ! IF ( ok_herbivores ) THEN d_ind(:,j) = d_ind(:,j) * EXP( - tau_eatup/herbivores(:,j) ) ENDIF ! ! 4.4 be sure that ind*cn_ind does not exceed 1 ! WHERE ( ( d_ind(:,j) .GT. zero ) .AND. & ( (ind(:,j)+d_ind(:,j))*cn_ind(:,j) .GT. un ) ) d_ind(:,j) = MAX( 1._r_std / cn_ind(:,j) - ind(:,j), 0._r_std ) ENDWHERE ! ! 4.5 new properties where there is establishment ( d_ind > 0 ) ! ! 4.5.1 biomass. ! Add biomass only if d_ind, over one year, is of the order of ind. ! (If we don't do this, the biomass density can become very low). ! In that case, take biomass from the atmosphere. ! compare establishment rate and present number of inidivuals many_new(:) = ( d_ind(:,j) .GT. 0.1 * ind(:,j) ) ! gives a better vectorization of the VPP IF ( ANY( many_new(:) ) ) THEN DO k = 1, nparts WHERE ( many_new(:) ) bm_new(:) = d_ind(:,j) * bm_sapl(j,k) / veget_max (:,j) biomass(:,j,k) = biomass(:,j,k) + bm_new(:) !NV passage 2D co2_to_bm(:,j) = co2_to_bm(:,j) + bm_new(:) / dt ENDWHERE ENDDO ! reset leaf ages. Should do a real calculation like in the npp routine, ! but this case is rare and not worth messing around. WHERE ( many_new(:) ) leaf_age(:,j,1) = zero leaf_frac(:,j,1) = un ENDWHERE DO m = 2, nleafages WHERE ( many_new(:) ) leaf_age(:,j,m) = zero leaf_frac(:,j,m) = zero ENDWHERE ENDDO ENDIF ! establishment rate is large WHERE ( d_ind(:,j) .GT. zero ) ! 4.5.2 age decreases age(:,j) = age(:,j) * ind(:,j) / ( ind(:,j) + d_ind(:,j) ) ! 4.5.3 new number of individuals ind(:,j) = ind(:,j) + d_ind(:,j) ENDWHERE ! ! 4.6 eventually convert excess sapwood to heartwood ! IF ( tree(j) ) THEN sm2(:) = zero WHERE ( d_ind(:,j) .GT. zero ) ! ratio of above / total sap parts sm_at(:) = biomass(:,j,isapabove) / & ( biomass(:,j,isapabove) + biomass(:,j,isapbelow) ) ! woodmass of an individual woodmass(:) = & ( biomass(:,j,isapabove) + biomass(:,j,isapbelow) + & biomass(:,j,iheartabove) + biomass(:,j,iheartbelow) ) / ind(:,j) ! crown area (m**2) depends on stem diameter (pipe model) dia(:) = ( woodmass(:) / ( pipe_density * pi/4. * pipe_tune2 ) ) & ** ( un / ( 2. + pipe_tune3 ) ) b1(:) = pipe_k1 / ( sla(j) * pipe_density*pipe_tune2 * dia(:)**pipe_tune3 ) * & ind(:,j) sm2(:) = lm_lastyearmax(:,j) / b1(:) ENDWHERE WHERE ( ( d_ind(:,j) .GT. zero ) .AND. & ( biomass(:,j,isapabove) + biomass(:,j,isapbelow) ) .GT. sm2(:) ) biomass(:,j,iheartabove) = biomass(:,j,iheartabove) + & ( biomass(:,j,isapabove) - sm2(:) * sm_at(:) ) biomass(:,j,isapabove) = sm2(:) * sm_at(:) biomass(:,j,iheartbelow) = biomass(:,j,iheartbelow) + & ( biomass(:,j,isapbelow) - sm2(:) * (un - sm_at) ) biomass(:,j,isapbelow) = sm2(:) * (un - sm_at(:)) ENDWHERE ENDIF ! tree ENDIF ! natural ENDDO ! loop over pfts ! ! 5 history ! d_ind = d_ind / dt CALL histwrite (hist_id_stomate, 'IND_ESTAB', itime, d_ind, npts*nvm, horipft_index) IF (bavard.GE.4) WRITE(numout,*) 'Leaving establish' END SUBROUTINE establish END MODULE lpj_establish