source: branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/lpj_cover.f90 @ 6319

Last change on this file since 6319 was 4647, checked in by josefine.ghattas, 7 years ago

Ticket #392 : Change variable name veget_max into veget_cov_max to avoid confusion with sechiba variable veget_max. The output variable VEGET_MAX is also changed into VEGET_COV_MAX.

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 18.2 KB
Line 
1! =================================================================================================================================
2! MODULE       : lpj_cover
3!
4! CONTACT      : orchidee-help _at_ listes.ipsl.fr
5!
6! LICENCE      : IPSL (2006)
7!                This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8!
9!>\BRIEF        Recalculate vegetation cover and LAI
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!!        Smith, B., I. C. Prentice, et al. (2001), Representation of vegetation
20!!        dynamics in the modelling of terrestrial ecosystems: comparing two
21!!        contrasting approaches within European climate space,
22!!        Global Ecology and Biogeography, 10, 621-637.\n
23!!
24!! SVN :
25!! $HeadURL$
26!! $Date$
27!! $Revision$
28!! \n
29!_ ================================================================================================================================
30
31MODULE lpj_cover
32
33  ! modules used:
34
35  USE ioipsl_para
36  USE stomate_data
37  USE pft_parameters
38
39  IMPLICIT NONE
40
41  ! private & public routines
42
43  PRIVATE
44  PUBLIC cover
45
46CONTAINS
47
48!! ================================================================================================================================
49!! SUBROUTINE     : lpj_cover
50!!
51!>\BRIEF          Recalculate vegetation cover and LAI
52!!
53!!\n DESCRIPTION : Veget_cov_max is first renewed here according to newly calculated foliage biomass in this calculation step
54!! Then, litter, soil carbon, and biomass are also recalcuted with taking into account the changes in Veget_cov_max (i.e. delta_veg)
55!! Grid-scale fpc (foliage projected coverage) is calculated to obtain the shadede ground area by leaf's light capture
56!! Finally, grid-scale fpc is adjusted not to exceed 1.0
57!!
58!! RECENT CHANGE(S) : None
59!!
60!! MAIN OUTPUT VARIABLE(S) : ::lai (leaf area index, @tex $(m^2 m^{-2})$ @endtex),
61!! :: veget (fractional vegetation cover, unitless)
62!!
63!! REFERENCE(S)   : None
64!!
65!! FLOWCHART :
66!! \latexonly
67!!     \includegraphics[scale=0.5]{lpj_cover_flowchart.png}
68!! \endlatexonly
69!! \n
70!_ ================================================================================================================================
71
72  SUBROUTINE cover (npts, cn_ind, ind, biomass, &
73       veget_cov_max, veget_cov_max_old, lai, litter, carbon, turnover_daily, bm_to_litter, &
74       co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, gpp_daily)
75
76!! 0. Variable and parameter declaration
77
78    !! 0.1 Input variables
79
80    INTEGER(i_std), INTENT(in)                                  :: npts             !! Domain size (unitless) 
81    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                :: cn_ind           !! Crown area
82                                                                                    !! @tex $(m^2)$ @endtex per individual
83    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                :: ind              !! Number of individuals
84                                                                                    !! @tex $(m^{-2})$ @endtex
85    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                :: veget_cov_max_old!! "Maximal" coverage fraction of a PFT (LAI->
86                                                                                    !! infinity) on ground at beginning of time
87
88    !! 0.2 Output variables
89
90    !! 0.3 Modified variables
91
92    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)             :: lai                 !! Leaf area index OF AN INDIVIDUAL PLANT
93                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
94    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: litter    !! Metabolic and structural litter, above and
95                                                                                       !! below ground @tex $(gC m^{-2})$ @endtex
96    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout)            :: carbon         !! Carbon pool: active, slow, or passive @tex $(gC m^{-2})$ @endtex
97    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass        !! Biomass @tex $(gC m^{-2})$ @endtex
98    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: veget_cov_max  !! "Maximal" coverage fraction of a PFT (LAI->
99                                                                                       !! infinity) on ground (unitless)
100    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: turnover_daily !! Turnover rates (gC m^{-2} day^{-1})
101    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter   !! Conversion of biomass to litter
102                                                                                       !! @tex $(gC m^{-2} day^{-1})$ @endtex
103    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: co2_to_bm      !! biomass up take for establishment           
104    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: co2_fire       !! Carbon emitted to the atmosphere by fire(living
105                                                                                       !! and dead biomass)
106    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: resp_hetero    !! Heterotrophic respiration
107    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: resp_maint     !! Maintenance respiration
108    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: resp_growth    !! Growth respiration
109    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: gpp_daily      !! Daily gross primary productivity
110
111    !! 0.4 Local variables
112
113    INTEGER(i_std)                                              :: i,j,k,m             !! Index (unitless)
114    REAL(r_std), DIMENSION(npts,nlitt,nlevs,nelements)          :: dilu_lit            !! Dilution for carbon variables
115    REAL(r_std), DIMENSION(npts,ncarb)                          :: dilu_soil_carbon    !! Dilution for carbon variables
116    REAL(r_std), DIMENSION(npts,nparts,nelements)               :: dilu_bio            !! Dilution for carbon variables
117    REAL(r_std), DIMENSION(npts)                                :: dilu_TCarbon        !! Dilution for carbon variables
118    REAL(r_std), DIMENSION(npts,nparts,nelements)               :: dilu_turnover_daily !! Dilution for carbon variables
119    REAL(r_std), DIMENSION(npts,nparts,nelements)               :: dilu_bm_to_litter   !! Dilution for carbon variables
120    REAL(r_std), DIMENSION(npts)                                :: dilu_co2flux_new    !! Dilution for carbon variables
121    REAL(r_std), DIMENSION(npts)                                :: dilu_gpp_daily      !! Dilution for carbon variables
122    REAL(r_std), DIMENSION(npts)                                :: dilu_resp_growth    !! Dilution for carbon variables
123    REAL(r_std), DIMENSION(npts)                                :: dilu_resp_maint     !! Dilution for carbon variables
124    REAL(r_std), DIMENSION(npts)                                :: dilu_resp_hetero    !! Dilution for carbon variables
125    REAL(r_std), DIMENSION(npts)                                :: dilu_co2_to_bm      !! Dilution for carbon variables
126    REAL(r_std), DIMENSION(npts)                                :: dilu_co2_fire       !! Dilution for carbon variables
127    REAL(r_std), DIMENSION(npts,nvm)                            :: TCarbon             !! Total carbon
128    REAL(r_std), DIMENSION(npts,nvm)                            :: co2flux_new         !! NBP after re-calculation in order to conserve carbon
129    REAL(r_std), DIMENSION(npts,nvm)                            :: co2flux_old         !! NBP before re-calculation
130    REAL(r_std), DIMENSION(nvm)                                 :: delta_veg           !! Conversion factors (unitless)
131    REAL(r_std), DIMENSION(nvm)                                 :: reduct              !! Conversion factors (unitless)
132    REAL(r_std)                                                 :: delta_veg_sum       !! Conversion factors (unitless)
133    REAL(r_std)                                                 :: diff                !! Conversion factors (unitless)
134    REAL(r_std)                                                 :: sr                  !! Conversion factors (unitless)
135    REAL(r_std), DIMENSION(npts)                                :: frac_nat            !! Conversion factors (unitless)
136    REAL(r_std), DIMENSION(npts)                                :: sum_vegettree       !! Conversion factors (unitless)
137    REAL(r_std), DIMENSION(npts)                                :: sum_vegetgrass      !! Conversion factors (unitless)
138    REAL(r_std), DIMENSION(npts)                                :: sum_veget_natveg    !! Conversion factors (unitless)
139    REAL(r_std), DIMENSION(npts)                                :: vartmp              !! Temporary variable used to add history
140
141!_ ================================================================================================================================
142
143 !! 1. If the vegetation is dynamic, calculate new maximum vegetation cover for natural plants
144 
145    IF ( ok_dgvm ) THEN
146
147       !! 1.1  Calculate initial values of vegetation cover
148       frac_nat(:) = un
149       sum_veget_natveg(:) = zero
150       veget_cov_max(:,ibare_sechiba) = un
151
152       DO j = 2,nvm ! loop over PFTs
153
154          IF ( natural(j) ) THEN
155             
156             ! Summation of individual tree crown area to get total foliar projected coverage
157             veget_cov_max(:,j) = ind(:,j) * cn_ind(:,j)
158             sum_veget_natveg(:) = sum_veget_natveg(:) + veget_cov_max(:,j)
159
160          ELSE
161             
162             !fraction occupied by agriculture needs to be substracted for the DGVM
163             !this is used below to constrain veget for natural vegetation, see below
164             frac_nat(:) = frac_nat(:) - veget_cov_max(:,j)
165
166          ENDIF
167
168       ENDDO ! loop over PFTs
169
170       DO i = 1, npts ! loop over grid points
171         
172          ! Recalculation of vegetation projected coverage when ::frac_nat was below ::sum_veget_natveg
173          ! It means that non-natural vegetation will recover ::veget_cov_max as natural vegetation
174          IF (sum_veget_natveg(i) .GT. frac_nat(i) .AND. frac_nat(i) .GT. min_stomate) THEN
175
176             DO j = 2,nvm ! loop over PFTs
177                IF( natural(j) ) THEN
178                   veget_cov_max(i,j) =  veget_cov_max(i,j) * frac_nat(i) / sum_veget_natveg(i)
179                ENDIF
180             ENDDO ! loop over PFTs
181
182          ENDIF
183       ENDDO ! loop over grid points
184       
185       ! Renew veget_cov_max of bare soil as 0 to difference of veget_cov_max (ibare_sechiba)
186       ! to current veget_cov_max
187       DO j = 2,nvm ! loop over PFTs
188          veget_cov_max(:,ibare_sechiba) = veget_cov_max(:,ibare_sechiba) - veget_cov_max(:,j)
189       ENDDO ! loop over PFTs
190       veget_cov_max(:,ibare_sechiba) = MAX( veget_cov_max(:,ibare_sechiba), zero )
191
192       !! 1.2 Calculate carbon fluxes between PFTs to maintain mass balance
193       !! Assure carbon closure when veget_cov_max changes(delta_veg): if veget_cov_max of some PFTs decrease, we use "dilu" to
194       !! record the corresponding lost in carbon (biomass, litter, soil carbon, gpp, respiration etc.) for
195       !! these PFTs, and re-allocate "dilu" to those PFTs with increasing veget_cov_max.
196       DO i = 1, npts ! loop over grid points
197         
198          ! Calculate the change in veget_cov_max between previous time step and current time step
199          delta_veg(:) = veget_cov_max(i,:)-veget_cov_max_old(i,:)
200          delta_veg_sum = SUM(delta_veg,MASK=delta_veg.LT.zero)
201
202          dilu_lit(i,:,:,:) = zero
203          dilu_soil_carbon(i,:) = zero
204          dilu_bio(i,:,:) = zero
205          dilu_TCarbon(i)=zero
206          dilu_turnover_daily(i,:,:)=zero
207          dilu_bm_to_litter(i,:,:)=zero
208          dilu_co2flux_new(i)=zero
209          dilu_gpp_daily(i)=zero
210          dilu_resp_growth(i)=zero
211          dilu_resp_maint(i)=zero
212          dilu_resp_hetero(i)=zero
213          dilu_co2_to_bm(i)=zero
214          dilu_co2_fire(i)=zero
215
216          ! Calculate TCarbon: total carbon including biomass, litter and soil carbon, as well as "today's" turnover and
217          ! bm_to_litter due to mortality, because today's turnover and bm_to_litter are not yet added into "litter" until tomorrow.
218          DO j=1, nvm
219                TCarbon(i,j)=SUM(biomass(i,j,:,:))+SUM(carbon(i,:,j))+SUM(litter(i,:,j,:,:))+&
220                     SUM(turnover_daily(i,j,:,:))+SUM(bm_to_litter(i,j,:,:))
221                co2flux_old(i,j)=resp_maint(i,j)+resp_growth(i,j)+resp_hetero(i,j)+co2_fire(i,j)-co2_to_bm(i,j)-gpp_daily(i,j)
222                co2flux_new(i,j)=resp_maint(i,j)+resp_growth(i,j)+resp_hetero(i,j)+co2_fire(i,j)-co2_to_bm(i,j)-gpp_daily(i,j)
223          ENDDO
224
225          DO j=1, nvm ! loop over PFTs
226             IF ( delta_veg(j) < -min_stomate ) THEN
227                dilu_lit(i,:,:,:) =  dilu_lit(i,:,:,:) + delta_veg(j) * litter(i,:,j,:,:) / delta_veg_sum
228                dilu_soil_carbon(i,:) =  dilu_soil_carbon(i,:) + delta_veg(j) * carbon(i,:,j) / delta_veg_sum
229                dilu_TCarbon(i)= dilu_TCarbon(i) + delta_veg(j) * TCarbon(i,j) / delta_veg_sum
230                dilu_turnover_daily(i,:,:)=dilu_turnover_daily(i,:,:)+delta_veg(j)*turnover_daily(i,j,:,:)/delta_veg_sum
231                dilu_bm_to_litter(i,:,:)=dilu_bm_to_litter(i,:,:)+delta_veg(j)*bm_to_litter(i,j,:,:)/delta_veg_sum
232                dilu_co2flux_new(i)=dilu_co2flux_new(i)+delta_veg(j)*co2flux_old(i,j)/delta_veg_sum
233                dilu_gpp_daily(i)=dilu_gpp_daily(i)+delta_veg(j)*gpp_daily(i,j)/delta_veg_sum
234                dilu_resp_growth(i)=dilu_resp_growth(i)+delta_veg(j)*resp_growth(i,j)/delta_veg_sum
235                dilu_resp_maint(i)=dilu_resp_maint(i)+delta_veg(j)*resp_maint(i,j)/delta_veg_sum
236                dilu_resp_hetero(i)=dilu_resp_hetero(i)+delta_veg(j)*resp_hetero(i,j)/delta_veg_sum
237                dilu_co2_to_bm(i)=dilu_co2_to_bm(i)+delta_veg(j)*co2_to_bm(i,j)/delta_veg_sum
238                dilu_co2_fire(i)=dilu_co2_fire(i)+delta_veg(j)*co2_fire(i,j)/delta_veg_sum
239             ENDIF
240          ENDDO ! loop over PFTs
241
242          DO j=1, nvm ! loop over PFTs
243             IF ( delta_veg(j) > min_stomate) THEN
244
245                ! Dilution of reservoirs
246                ! Recalculate the litter and soil carbon with taking into accout the change in
247                ! veget_cov_max (delta_veg)
248                ! Litter
249                litter(i,:,j,:,:)=(litter(i,:,j,:,:) * veget_cov_max_old(i,j) + dilu_lit(i,:,:,:) * delta_veg(j)) &
250                                  / veget_cov_max(i,j)
251
252                ! Soil carbon
253                carbon(i,:,j)=(carbon(i,:,j) * veget_cov_max_old(i,j) + dilu_soil_carbon(i,:) * delta_veg(j)) / veget_cov_max(i,j)
254                TCarbon(i,j)=(TCarbon(i,j) * veget_cov_max_old(i,j) + dilu_TCarbon(i) * delta_veg(j)) / veget_cov_max(i,j)
255                turnover_daily(i,j,:,:)=(turnover_daily(i,j,:,:)*veget_cov_max_old(i,j)+&
256                     dilu_turnover_daily(i,:,:)*delta_veg(j))/veget_cov_max(i,j)
257                bm_to_litter(i,j,:,:)=(bm_to_litter(i,j,:,:)*veget_cov_max_old(i,j)+&
258                     dilu_bm_to_litter(i,:,:)*delta_veg(j))/veget_cov_max(i,j)
259                co2flux_new(i,j)=(co2flux_old(i,j)*veget_cov_max_old(i,j)+dilu_co2flux_new(i)*delta_veg(j))/veget_cov_max(i,j)
260                gpp_daily(i,j)=(gpp_daily(i,j)*veget_cov_max_old(i,j)+dilu_gpp_daily(i)*delta_veg(j))/veget_cov_max(i,j)
261                resp_growth(i,j)=(resp_growth(i,j)*veget_cov_max_old(i,j)+dilu_resp_growth(i)*delta_veg(j))/veget_cov_max(i,j)
262                resp_maint(i,j)=(resp_maint(i,j)*veget_cov_max_old(i,j)+dilu_resp_maint(i)*delta_veg(j))/veget_cov_max(i,j)
263                resp_hetero(i,j)=(resp_hetero(i,j)*veget_cov_max_old(i,j)+dilu_resp_hetero(i)*delta_veg(j))/veget_cov_max(i,j)
264                co2_to_bm(i,j)=(co2_to_bm(i,j)*veget_cov_max_old(i,j)+dilu_co2_to_bm(i)*delta_veg(j))/veget_cov_max(i,j)
265                co2_fire(i,j)=(co2_fire(i,j)*veget_cov_max_old(i,j)+dilu_co2_fire(i)*delta_veg(j))/veget_cov_max(i,j)
266             ENDIF
267
268             IF(veget_cov_max(i,j).GT.min_stomate) THEN
269
270                ! Correct biomass densities to conserve mass
271                ! since it's defined on veget_cov_max
272                biomass(i,j,:,:) = biomass(i,j,:,:) * veget_cov_max_old(i,j) / veget_cov_max(i,j)
273
274             ENDIF
275
276          ENDDO ! loop over PFTs
277       ENDDO ! loop over grid points
278
279       vartmp(:)=SUM(gpp_daily*veget_cov_max,dim=2)
280       CALL histwrite_p (hist_id_stomate, "tGPP", itime, vartmp, npts, hori_index)
281       vartmp(:)=SUM(resp_growth*veget_cov_max,dim=2)
282       CALL histwrite_p (hist_id_stomate, "tRESP_GROWTH", itime, vartmp, npts, hori_index)
283       vartmp(:)=SUM(resp_maint*veget_cov_max,dim=2)
284       CALL histwrite_p (hist_id_stomate, "tRESP_MAINT", itime, vartmp, npts, hori_index)
285       vartmp(:)=SUM(resp_hetero*veget_cov_max,dim=2)
286       CALL histwrite_p (hist_id_stomate, "tRESP_HETERO", itime, vartmp, npts, hori_index)
287       vartmp(:)=SUM(co2_to_bm*veget_cov_max,dim=2)
288       CALL histwrite_p (hist_id_stomate, "tCO2_TAKEN", itime, vartmp, npts, hori_index)
289       vartmp(:)=SUM(co2_fire*veget_cov_max,dim=2)
290       CALL histwrite_p (hist_id_stomate, "tCO2_FIRE", itime, vartmp, npts, hori_index)
291       vartmp(:)=SUM(co2flux_new*veget_cov_max,dim=2)
292       CALL histwrite_p (hist_id_stomate, "tCO2FLUX", itime, vartmp, npts, hori_index)
293       vartmp(:)=SUM(co2flux_old*veget_cov_max_old,dim=2)
294       CALL histwrite_p (hist_id_stomate, "tCO2FLUX_OLD", itime, vartmp, npts, hori_index)
295       vartmp(:)=SUM(TCarbon*veget_cov_max,dim=2)
296       CALL histwrite_p (hist_id_stomate, "tCARBON", itime, vartmp, npts, hori_index)
297       vartmp(:)=SUM(SUM(biomass(:,:,:,icarbon),dim=3)*veget_cov_max,dim=2)
298       CALL histwrite_p (hist_id_stomate, "tBIOMASS", itime, vartmp, npts, hori_index)
299       vartmp(:)=SUM(SUM(SUM(litter(:,:,:,:,icarbon),dim=4),dim=2)*veget_cov_max,dim=2)
300       CALL histwrite_p (hist_id_stomate, "tLITTER", itime, vartmp, npts, hori_index)
301       vartmp(:)=SUM(SUM(carbon,dim=2)*veget_cov_max,dim=2)
302       CALL histwrite_p (hist_id_stomate, "tSOILC", itime, vartmp, npts, hori_index)
303    ENDIF
304
305  END SUBROUTINE cover
306
307END MODULE lpj_cover
Note: See TracBrowser for help on using the repository browser.