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

Last change on this file since 6393 was 6369, checked in by josefine.ghattas, 5 years ago

Integrated correction of diagnostic variables rhSoil and rhLitter as done in the trunk rev [6362]. See ticket #630

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 19.3 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_hetero_litter, resp_hetero_soil, 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_hetero_litter  !! Heterotrophic respiration from litter
108    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: resp_hetero_soil    !! Heterotrophic respiration from soil
109    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: resp_maint     !! Maintenance respiration
110    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: resp_growth    !! Growth respiration
111    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                  :: gpp_daily      !! Daily gross primary productivity
112
113    !! 0.4 Local variables
114
115    INTEGER(i_std)                                              :: i,j,k,m             !! Index (unitless)
116    REAL(r_std), DIMENSION(npts,nlitt,nlevs,nelements)          :: dilu_lit            !! Dilution for carbon variables
117    REAL(r_std), DIMENSION(npts,ncarb)                          :: dilu_soil_carbon    !! Dilution for carbon variables
118    REAL(r_std), DIMENSION(npts,nparts,nelements)               :: dilu_bio            !! Dilution for carbon variables
119    REAL(r_std), DIMENSION(npts)                                :: dilu_TCarbon        !! Dilution for carbon variables
120    REAL(r_std), DIMENSION(npts,nparts,nelements)               :: dilu_turnover_daily !! Dilution for carbon variables
121    REAL(r_std), DIMENSION(npts,nparts,nelements)               :: dilu_bm_to_litter   !! Dilution for carbon variables
122    REAL(r_std), DIMENSION(npts)                                :: dilu_co2flux_new    !! Dilution for carbon variables
123    REAL(r_std), DIMENSION(npts)                                :: dilu_gpp_daily      !! Dilution for carbon variables
124    REAL(r_std), DIMENSION(npts)                                :: dilu_resp_growth    !! Dilution for carbon variables
125    REAL(r_std), DIMENSION(npts)                                :: dilu_resp_maint     !! Dilution for carbon variables
126    REAL(r_std), DIMENSION(npts)                                :: dilu_resp_hetero    !! Dilution for carbon variables
127    REAL(r_std), DIMENSION(npts)                                :: dilu_resp_hetero_litter  !! Dilution for carbon variables
128    REAL(r_std), DIMENSION(npts)                                :: dilu_resp_hetero_soil    !! Dilution for carbon variables
129    REAL(r_std), DIMENSION(npts)                                :: dilu_co2_to_bm      !! Dilution for carbon variables
130    REAL(r_std), DIMENSION(npts)                                :: dilu_co2_fire       !! Dilution for carbon variables
131    REAL(r_std), DIMENSION(npts,nvm)                            :: TCarbon             !! Total carbon
132    REAL(r_std), DIMENSION(npts,nvm)                            :: co2flux_new         !! NBP after re-calculation in order to conserve carbon
133    REAL(r_std), DIMENSION(npts,nvm)                            :: co2flux_old         !! NBP before re-calculation
134    REAL(r_std), DIMENSION(nvm)                                 :: delta_veg           !! Conversion factors (unitless)
135    REAL(r_std), DIMENSION(nvm)                                 :: reduct              !! Conversion factors (unitless)
136    REAL(r_std)                                                 :: delta_veg_sum       !! Conversion factors (unitless)
137    REAL(r_std)                                                 :: diff                !! Conversion factors (unitless)
138    REAL(r_std)                                                 :: sr                  !! Conversion factors (unitless)
139    REAL(r_std), DIMENSION(npts)                                :: frac_nat            !! Conversion factors (unitless)
140    REAL(r_std), DIMENSION(npts)                                :: sum_vegettree       !! Conversion factors (unitless)
141    REAL(r_std), DIMENSION(npts)                                :: sum_vegetgrass      !! Conversion factors (unitless)
142    REAL(r_std), DIMENSION(npts)                                :: sum_veget_natveg    !! Conversion factors (unitless)
143    REAL(r_std), DIMENSION(npts)                                :: vartmp              !! Temporary variable used to add history
144
145!_ ================================================================================================================================
146
147 !! 1. If the vegetation is dynamic, calculate new maximum vegetation cover for natural plants
148 
149    IF ( ok_dgvm ) THEN
150
151       !! 1.1  Calculate initial values of vegetation cover
152       frac_nat(:) = un
153       sum_veget_natveg(:) = zero
154       veget_cov_max(:,ibare_sechiba) = un
155
156       DO j = 2,nvm ! loop over PFTs
157
158          IF ( natural(j) ) THEN
159             
160             ! Summation of individual tree crown area to get total foliar projected coverage
161             veget_cov_max(:,j) = ind(:,j) * cn_ind(:,j)
162             sum_veget_natveg(:) = sum_veget_natveg(:) + veget_cov_max(:,j)
163
164          ELSE
165             
166             !fraction occupied by agriculture needs to be substracted for the DGVM
167             !this is used below to constrain veget for natural vegetation, see below
168             frac_nat(:) = frac_nat(:) - veget_cov_max(:,j)
169
170          ENDIF
171
172       ENDDO ! loop over PFTs
173
174       DO i = 1, npts ! loop over grid points
175         
176          ! Recalculation of vegetation projected coverage when ::frac_nat was below ::sum_veget_natveg
177          ! It means that non-natural vegetation will recover ::veget_cov_max as natural vegetation
178          IF (sum_veget_natveg(i) .GT. frac_nat(i) .AND. frac_nat(i) .GT. min_stomate) THEN
179
180             DO j = 2,nvm ! loop over PFTs
181                IF( natural(j) ) THEN
182                   veget_cov_max(i,j) =  veget_cov_max(i,j) * frac_nat(i) / sum_veget_natveg(i)
183                ENDIF
184             ENDDO ! loop over PFTs
185
186          ENDIF
187       ENDDO ! loop over grid points
188       
189       ! Renew veget_cov_max of bare soil as 0 to difference of veget_cov_max (ibare_sechiba)
190       ! to current veget_cov_max
191       DO j = 2,nvm ! loop over PFTs
192          veget_cov_max(:,ibare_sechiba) = veget_cov_max(:,ibare_sechiba) - veget_cov_max(:,j)
193       ENDDO ! loop over PFTs
194       veget_cov_max(:,ibare_sechiba) = MAX( veget_cov_max(:,ibare_sechiba), zero )
195
196       !! 1.2 Calculate carbon fluxes between PFTs to maintain mass balance
197       !! Assure carbon closure when veget_cov_max changes(delta_veg): if veget_cov_max of some PFTs decrease, we use "dilu" to
198       !! record the corresponding lost in carbon (biomass, litter, soil carbon, gpp, respiration etc.) for
199       !! these PFTs, and re-allocate "dilu" to those PFTs with increasing veget_cov_max.
200       DO i = 1, npts ! loop over grid points
201         
202          ! Calculate the change in veget_cov_max between previous time step and current time step
203          delta_veg(:) = veget_cov_max(i,:)-veget_cov_max_old(i,:)
204          delta_veg_sum = SUM(delta_veg,MASK=delta_veg.LT.zero)
205
206          dilu_lit(i,:,:,:) = zero
207          dilu_soil_carbon(i,:) = zero
208          dilu_bio(i,:,:) = zero
209          dilu_TCarbon(i)=zero
210          dilu_turnover_daily(i,:,:)=zero
211          dilu_bm_to_litter(i,:,:)=zero
212          dilu_co2flux_new(i)=zero
213          dilu_gpp_daily(i)=zero
214          dilu_resp_growth(i)=zero
215          dilu_resp_maint(i)=zero
216          dilu_resp_hetero(i)=zero
217          dilu_resp_hetero_litter(i)=zero
218          dilu_resp_hetero_soil(i)=zero
219          dilu_co2_to_bm(i)=zero
220          dilu_co2_fire(i)=zero
221
222          ! Calculate TCarbon: total carbon including biomass, litter and soil carbon, as well as "today's" turnover and
223          ! bm_to_litter due to mortality, because today's turnover and bm_to_litter are not yet added into "litter" until tomorrow.
224          DO j=1, nvm
225                TCarbon(i,j)=SUM(biomass(i,j,:,:))+SUM(carbon(i,:,j))+SUM(litter(i,:,j,:,:))+&
226                     SUM(turnover_daily(i,j,:,:))+SUM(bm_to_litter(i,j,:,:))
227                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)
228                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)
229          ENDDO
230
231          DO j=1, nvm ! loop over PFTs
232             IF ( delta_veg(j) < -min_stomate ) THEN
233                dilu_lit(i,:,:,:) =  dilu_lit(i,:,:,:) + delta_veg(j) * litter(i,:,j,:,:) / delta_veg_sum
234                dilu_soil_carbon(i,:) =  dilu_soil_carbon(i,:) + delta_veg(j) * carbon(i,:,j) / delta_veg_sum
235                dilu_TCarbon(i)= dilu_TCarbon(i) + delta_veg(j) * TCarbon(i,j) / delta_veg_sum
236                dilu_turnover_daily(i,:,:)=dilu_turnover_daily(i,:,:)+delta_veg(j)*turnover_daily(i,j,:,:)/delta_veg_sum
237                dilu_bm_to_litter(i,:,:)=dilu_bm_to_litter(i,:,:)+delta_veg(j)*bm_to_litter(i,j,:,:)/delta_veg_sum
238                dilu_co2flux_new(i)=dilu_co2flux_new(i)+delta_veg(j)*co2flux_old(i,j)/delta_veg_sum
239                dilu_gpp_daily(i)=dilu_gpp_daily(i)+delta_veg(j)*gpp_daily(i,j)/delta_veg_sum
240                dilu_resp_growth(i)=dilu_resp_growth(i)+delta_veg(j)*resp_growth(i,j)/delta_veg_sum
241                dilu_resp_maint(i)=dilu_resp_maint(i)+delta_veg(j)*resp_maint(i,j)/delta_veg_sum
242                dilu_resp_hetero(i)=dilu_resp_hetero(i)+delta_veg(j)*resp_hetero(i,j)/delta_veg_sum
243                dilu_resp_hetero_litter(i)=dilu_resp_hetero_litter(i)+delta_veg(j)*resp_hetero_litter(i,j)/delta_veg_sum
244                dilu_resp_hetero_soil(i)=dilu_resp_hetero_soil(i)+delta_veg(j)*resp_hetero_soil(i,j)/delta_veg_sum
245                dilu_co2_to_bm(i)=dilu_co2_to_bm(i)+delta_veg(j)*co2_to_bm(i,j)/delta_veg_sum
246                dilu_co2_fire(i)=dilu_co2_fire(i)+delta_veg(j)*co2_fire(i,j)/delta_veg_sum
247             ENDIF
248          ENDDO ! loop over PFTs
249
250          DO j=1, nvm ! loop over PFTs
251             IF ( delta_veg(j) > min_stomate) THEN
252
253                ! Dilution of reservoirs
254                ! Recalculate the litter and soil carbon with taking into accout the change in
255                ! veget_cov_max (delta_veg)
256                ! Litter
257                litter(i,:,j,:,:)=(litter(i,:,j,:,:) * veget_cov_max_old(i,j) + dilu_lit(i,:,:,:) * delta_veg(j)) &
258                                  / veget_cov_max(i,j)
259
260                ! Soil carbon
261                carbon(i,:,j)=(carbon(i,:,j) * veget_cov_max_old(i,j) + dilu_soil_carbon(i,:) * delta_veg(j)) / veget_cov_max(i,j)
262                TCarbon(i,j)=(TCarbon(i,j) * veget_cov_max_old(i,j) + dilu_TCarbon(i) * delta_veg(j)) / veget_cov_max(i,j)
263                turnover_daily(i,j,:,:)=(turnover_daily(i,j,:,:)*veget_cov_max_old(i,j)+&
264                     dilu_turnover_daily(i,:,:)*delta_veg(j))/veget_cov_max(i,j)
265                bm_to_litter(i,j,:,:)=(bm_to_litter(i,j,:,:)*veget_cov_max_old(i,j)+&
266                     dilu_bm_to_litter(i,:,:)*delta_veg(j))/veget_cov_max(i,j)
267                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)
268                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)
269                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)
270                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)
271                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)
272                resp_hetero_litter(i,j)=(resp_hetero_litter(i,j)*veget_cov_max_old(i,j)+ &
273                       dilu_resp_hetero_litter(i)*delta_veg(j))/veget_cov_max(i,j)
274                resp_hetero_soil(i,j)=(resp_hetero_soil(i,j)*veget_cov_max_old(i,j)+ &
275                       dilu_resp_hetero_soil(i)*delta_veg(j))/veget_cov_max(i,j)
276                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)
277                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)
278             ENDIF
279
280             IF(veget_cov_max(i,j).GT.min_stomate) THEN
281
282                ! Correct biomass densities to conserve mass
283                ! since it's defined on veget_cov_max
284                biomass(i,j,:,:) = biomass(i,j,:,:) * veget_cov_max_old(i,j) / veget_cov_max(i,j)
285
286             ENDIF
287
288          ENDDO ! loop over PFTs
289       ENDDO ! loop over grid points
290
291       vartmp(:)=SUM(gpp_daily*veget_cov_max,dim=2)
292       CALL histwrite_p (hist_id_stomate, "tGPP", itime, vartmp, npts, hori_index)
293       vartmp(:)=SUM(resp_growth*veget_cov_max,dim=2)
294       CALL histwrite_p (hist_id_stomate, "tRESP_GROWTH", itime, vartmp, npts, hori_index)
295       vartmp(:)=SUM(resp_maint*veget_cov_max,dim=2)
296       CALL histwrite_p (hist_id_stomate, "tRESP_MAINT", itime, vartmp, npts, hori_index)
297       vartmp(:)=SUM(resp_hetero*veget_cov_max,dim=2)
298       CALL histwrite_p (hist_id_stomate, "tRESP_HETERO", itime, vartmp, npts, hori_index)
299       vartmp(:)=SUM(co2_to_bm*veget_cov_max,dim=2)
300       CALL histwrite_p (hist_id_stomate, "tCO2_TAKEN", itime, vartmp, npts, hori_index)
301       vartmp(:)=SUM(co2_fire*veget_cov_max,dim=2)
302       CALL histwrite_p (hist_id_stomate, "tCO2_FIRE", itime, vartmp, npts, hori_index)
303       vartmp(:)=SUM(co2flux_new*veget_cov_max,dim=2)
304       CALL histwrite_p (hist_id_stomate, "tCO2FLUX", itime, vartmp, npts, hori_index)
305       vartmp(:)=SUM(co2flux_old*veget_cov_max_old,dim=2)
306       CALL histwrite_p (hist_id_stomate, "tCO2FLUX_OLD", itime, vartmp, npts, hori_index)
307       vartmp(:)=SUM(TCarbon*veget_cov_max,dim=2)
308       CALL histwrite_p (hist_id_stomate, "tCARBON", itime, vartmp, npts, hori_index)
309       vartmp(:)=SUM(SUM(biomass(:,:,:,icarbon),dim=3)*veget_cov_max,dim=2)
310       CALL histwrite_p (hist_id_stomate, "tBIOMASS", itime, vartmp, npts, hori_index)
311       vartmp(:)=SUM(SUM(SUM(litter(:,:,:,:,icarbon),dim=4),dim=2)*veget_cov_max,dim=2)
312       CALL histwrite_p (hist_id_stomate, "tLITTER", itime, vartmp, npts, hori_index)
313       vartmp(:)=SUM(SUM(carbon,dim=2)*veget_cov_max,dim=2)
314       CALL histwrite_p (hist_id_stomate, "tSOILC", itime, vartmp, npts, hori_index)
315    ENDIF
316
317  END SUBROUTINE cover
318
319END MODULE lpj_cover
Note: See TracBrowser for help on using the repository browser.