source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_stomate/lpj_cover.f90 @ 5816

Last change on this file since 5816 was 5816, checked in by jinfeng.chang, 5 years ago

copy ORCHIDEE-GMv3.2 for publication

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