Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lpj_cover.f90 @ 8624

Last change on this file since 8624 was 3549, checked in by bertrand.guenet, 8 years ago
update between the trunk rev 3340 and SOM it's still not working
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 20.3 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
31	MODULE lpj_cover
32
33	! modules used:
34
35	USE ioipsl_para
36	USE stomate_data
37	USE pft_parameters
38	USE constantes_soil
39
40	IMPLICIT NONE
41
42	! private & public routines
43
44	PRIVATE
45	PUBLIC cover
46
47	CONTAINS
48
49	!! ================================================================================================================================
50	!! SUBROUTINE : lpj_cover
51	!!
52	!>\BRIEF Recalculate vegetation cover and LAI
53	!!
54	!!\n DESCRIPTION : Veget_max is first renewed here according to newly calculated foliage biomass in this calculation step
55	!! Then, litter, soil carbon, and biomass are also recalcuted with taking into account the changes in Veget_max (i.e. delta_veg)
56	!! Grid-scale fpc (foliage projected coverage) is calculated to obtain the shadede ground area by leaf's light capture
57	!! Finally, grid-scale fpc is adjusted not to exceed 1.0
58	!!
59	!! RECENT CHANGE(S) : None
60	!!
61	!! MAIN OUTPUT VARIABLE(S) : ::lai (leaf area index, @tex $(m^2 m^{-2})$ @endtex),
62	!! :: veget (fractional vegetation cover, unitless)
63	!!
64	!! REFERENCE(S) : None
65	!!
66	!! FLOWCHART :
67	!! \latexonly
68	!! \includegraphics[scale=0.5]{lpj_cover_flowchart.png}
69	!! \endlatexonly
70	!! \n
71	!_ ================================================================================================================================
72
73	SUBROUTINE cover (npts, cn_ind, ind, biomass, &
74	veget_max, veget_max_old, lai, litter_above, litter_below, carbon, DOC, turnover_daily, bm_to_litter, &
75	co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, gpp_daily, moist_soil)
76
77	!! 0. Variable and parameter declaration
78
79	!! 0.1 Input variables
80
81	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
82	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: cn_ind !! Crown area
83	!! @tex $(m^2)$ @endtex per individual
84	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: ind !! Number of individuals
85	!! @tex $(m^{-2})$ @endtex
86	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max_old !! "Maximal" coverage fraction of a PFT (LAI->
87	!! infinity) on ground at beginning of time
88	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: moist_soil !! moisture in soil for one pixel (m3/m3)
89	!! 0.2 Output variables
90
91	!! 0.3 Modified variables
92
93	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lai !! Leaf area index OF AN INDIVIDUAL PLANT
94	!! @tex $(m^2 m^{-2})$ @endtex
95	REAL(r_std), DIMENSION(npts,nlitt,nvm,nelements), INTENT(inout) :: litter_above !! Metabolic and structural litter, above
96	!! ground @tex $(gC m^{-2})$ @endtex
97	REAL(r_std), DIMENSION(npts,nlitt,nvm,nbdl,nelements), INTENT(inout) :: litter_below !! Metabolic and structural litter, below ground
98	!! @tex $(gC m^{-2})$ @endtex
99	REAL(r_std), DIMENSION(npts,ncarb,nvm,nbdl), INTENT(inout) :: carbon !! Carbon pool: active, slow, or passive
100	!! @tex $(gC m^{-2})$ @endtex
101	REAL(r_std), DIMENSION(npts,nvm,nbdl,ndoc,npool,nelements), INTENT(inout) :: DOC !! Dissolved Organic Carbon in soil
102	!! The unit is given by m^3 of water
103	!! @tex $(gC m{-3} of water) $ @endtex
104	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass !! Biomass @tex $(gC m^{-2})$ @endtex
105	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max !! "Maximal" coverage fraction of a PFT (LAI->
106	!! infinity) on ground (unitless)
107	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: turnover_daily !! Turnover rates (gC m^{-2} day^{-1})
108	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! Conversion of biomass to litter
109	!! @tex $(gC m^{-2} day^{-1})$ @endtex
110	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_to_bm !! biomass up take for establishment
111	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_fire !! Carbon emitted to the atmosphere by fire(living
112	!! and dead biomass)
113	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
114	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_maint !! Maintenance respiration
115	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_growth !! Growth respiration
116	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
117
118	!! 0.4 Local variables
119
120	INTEGER(i_std) :: i,j,k,m !! Index (unitless)
121	REAL(r_std), DIMENSION(npts,nlitt,nelements) :: dilu_lit_above !! Above ground litter dilution @tex $(gC m^{-2})$ @endtex
122	REAL(r_std), DIMENSION(npts,nlitt,nbdl,nelements) :: dilu_lit_below !! Below ground litter dilution @tex $(gC m^{-2})$ @endtex
123	REAL(r_std), DIMENSION(npts,ncarb,nbdl) :: dilu_soil_carbon !! Dilution for carbon variables
124	REAL(r_std), DIMENSION(npts,nparts,nelements) :: dilu_bio !! Dilution for carbon variables
125	REAL(r_std), DIMENSION(npts) :: dilu_TCarbon !! Dilution for carbon variables
126	REAL(r_std), DIMENSION(npts,nparts,nelements) :: dilu_turnover_daily !! Dilution for carbon variables
127	REAL(r_std), DIMENSION(npts,nparts,nelements) :: dilu_bm_to_litter !! Dilution for carbon variables
128	REAL(r_std), DIMENSION(npts) :: dilu_co2flux_new !! Dilution for carbon variables
129	REAL(r_std), DIMENSION(npts) :: dilu_gpp_daily !! Dilution for carbon variables
130	REAL(r_std), DIMENSION(npts) :: dilu_resp_growth !! Dilution for carbon variables
131	REAL(r_std), DIMENSION(npts) :: dilu_resp_maint !! Dilution for carbon variables
132	REAL(r_std), DIMENSION(npts) :: dilu_resp_hetero !! Dilution for carbon variables
133	REAL(r_std), DIMENSION(npts) :: dilu_co2_to_bm !! Dilution for carbon variables
134	REAL(r_std), DIMENSION(npts) :: dilu_co2_fire !! Dilution for carbon variables
135	REAL(r_std), DIMENSION(npts,nvm) :: TCarbon !! Total carbon
136	REAL(r_std), DIMENSION(npts,nvm) :: co2flux_new !! NBP after re-calculation in order to conserve carbon
137	REAL(r_std), DIMENSION(npts,nvm) :: co2flux_old !! NBP before re-calculation
138	REAL(r_std), DIMENSION(nvm) :: delta_veg !! Conversion factors (unitless)
139	REAL(r_std), DIMENSION(nvm) :: reduct !! Conversion factors (unitless)
140	REAL(r_std) :: delta_veg_sum !! Conversion factors (unitless)
141	REAL(r_std) :: diff !! Conversion factors (unitless)
142	REAL(r_std) :: sr !! Conversion factors (unitless)
143	REAL(r_std), DIMENSION(npts) :: frac_nat !! Conversion factors (unitless)
144	REAL(r_std), DIMENSION(npts) :: sum_vegettree !! Conversion factors (unitless)
145	REAL(r_std), DIMENSION(npts) :: sum_vegetgrass !! Conversion factors (unitless)
146	REAL(r_std), DIMENSION(npts) :: sum_veget_natveg !! Conversion factors (unitless)
147	REAL(r_std), DIMENSION(npts) :: vartmp !! Temporary variable used to add history
148	REAL(r_std), DIMENSION(0:nbdl) :: z_soil !! Variable to store depth of the different
149	!! soil layers (m)
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
200	veget_max(:,ibare_sechiba) = MAX( veget_max(:,ibare_sechiba), zero )
201	!! 1.1.1 Copy the depth of the different soil layers (number of layers=nbdl)
202	! previously calculated as variable diaglev in routines
203	! sechiba.f90 and slowproc.f90
204	z_soil(0) = zero
205	z_soil(1:nbdl) = diaglev(1:nbdl)
206
207	!! 1.2 Calculate carbon fluxes between PFTs to maintain mass balance
208	!! Assure carbon closure when veget_max changes(delta_veg): if veget_max of some PFTs decrease, we use "dilu" to
209	!! record the corresponding lost in carbon (biomass, litter, soil carbon, gpp, respiration etc.) for
210	!! these PFTs, and re-allocate "dilu" to those PFTs with increasing veget_max.
211	DO i = 1, npts ! loop over grid points
212
213	! Calculate the change in veget_max between previous time step and current time step
214	delta_veg(:) = veget_max(i,:)-veget_max_old(i,:)
215	delta_veg_sum = SUM(delta_veg,MASK=delta_veg.LT.zero)
216
217	dilu_lit_above(i,:,:) = zero
218	dilu_lit_below(i,:,:,:) = zero
219	dilu_soil_carbon(i,:,:) = zero
220
221	dilu_bio(i,:,:) = zero
222	dilu_TCarbon(i)=zero
223	dilu_turnover_daily(i,:,:)=zero
224	dilu_bm_to_litter(i,:,:)=zero
225	dilu_co2flux_new(i)=zero
226	dilu_gpp_daily(i)=zero
227	dilu_resp_growth(i)=zero
228	dilu_resp_maint(i)=zero
229	dilu_resp_hetero(i)=zero
230	dilu_co2_to_bm(i)=zero
231	dilu_co2_fire(i)=zero
232
233	! Calculate TCarbon: total carbon including biomass, litter and soil carbon, as well as "today's" turnover and
234	! bm_to_litter due to mortality, because today's turnover and bm_to_litter are not yet added into "litter" until tomorrow.
235	DO j=1, nvm
236	TCarbon(i,j)=SUM(biomass(i,j,:,:))+SUM(carbon(i,:,j,:))+SUM(litter_above(i,:,j,:))+ SUM(litter_below(i,:,j,:,:))+&
237	SUM(turnover_daily(i,j,:,:))+SUM(bm_to_litter(i,j,:,:))
238	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)
239	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)
240	ENDDO
241
242	DO j=1, nvm ! loop over PFTs
243	DO k=2,nbdl
244	IF ( delta_veg(j) < -min_stomate ) THEN
245	dilu_lit_above(i,:,:) = dilu_lit_above(i,:,:) + delta_veg(j) * litter_above(i,:,j,:) / delta_veg_sum
246	dilu_lit_below(i,:,:,:) = dilu_lit_below(i,:,:,:) + delta_veg(j) * litter_below(i,:,j,:,:) / delta_veg_sum
247	DO m = 1,npool
248	dilu_soil_carbon(i,:,k) = dilu_soil_carbon(i,:,k) + delta_veg(j) * (carbon(i,:,j,k)+ &
249	(DOC(:,j,k,ifree,m,icarbon) + DOC(:,j,k,iadsorbed,m,icarbon))* moist_soil(:,k) * (z_soil(k)-z_soil(k-1))) / delta_veg_sum
250	ENDDO
251	dilu_TCarbon(i)= dilu_TCarbon(i) + delta_veg(j) * TCarbon(i,j) / delta_veg_sum
252	dilu_turnover_daily(i,:,:)=dilu_turnover_daily(i,:,:)+delta_veg(j)*turnover_daily(i,j,:,:)/delta_veg_sum
253	dilu_bm_to_litter(i,:,:)=dilu_bm_to_litter(i,:,:)+delta_veg(j)*bm_to_litter(i,j,:,:)/delta_veg_sum
254	dilu_co2flux_new(i)=dilu_co2flux_new(i)+delta_veg(j)*co2flux_old(i,j)/delta_veg_sum
255	dilu_gpp_daily(i)=dilu_gpp_daily(i)+delta_veg(j)*gpp_daily(i,j)/delta_veg_sum
256	dilu_resp_growth(i)=dilu_resp_growth(i)+delta_veg(j)*resp_growth(i,j)/delta_veg_sum
257	dilu_resp_maint(i)=dilu_resp_maint(i)+delta_veg(j)*resp_maint(i,j)/delta_veg_sum
258	dilu_resp_hetero(i)=dilu_resp_hetero(i)+delta_veg(j)*resp_hetero(i,j)/delta_veg_sum
259	dilu_co2_to_bm(i)=dilu_co2_to_bm(i)+delta_veg(j)*co2_to_bm(i,j)/delta_veg_sum
260	dilu_co2_fire(i)=dilu_co2_fire(i)+delta_veg(j)*co2_fire(i,j)/delta_veg_sum
261	ENDIF
262	ENDDO
263	DO m =1,npool
264	dilu_soil_carbon(i,:,1) = dilu_soil_carbon(i,:,1) + delta_veg(j) * (carbon(i,:,j,1)+ &
265	(DOC(:,j,1,ifree,m,icarbon) + DOC(:,j,1,iadsorbed,m,icarbon))* moist_soil(:,1) * (z_soil(1))) / delta_veg_sum
266	ENDDO
267	ENDDO ! loop over PFTs
268
269	DO j=1, nvm ! loop over PFTs
270	IF ( delta_veg(j) > min_stomate) THEN
271
272	! Dilution of reservoirs
273	! Recalculate the litter and soil carbon with taking into accout the change in
274	! veget_max (delta_veg)
275	! Litter
276	litter_above(i,:,j,:)=(litter_above(i,:,j,:) * veget_max_old(i,j) + dilu_lit_above(i,:,:) * delta_veg(j)) / veget_max(i,j)
277	litter_below(i,:,j,:,:)=(litter_below(i,:,j,:,:) * veget_max_old(i,j) + dilu_lit_below(i,:,:,:) * delta_veg(j)) / veget_max(i,j)
278	! Soil carbon
279	carbon(i,:,j,:)=(carbon(i,:,j,:) * veget_max_old(i,j) + dilu_soil_carbon(i,:,:) * delta_veg(j)) / veget_max(i,j)
280	TCarbon(i,j)=(TCarbon(i,j) * veget_max_old(i,j) + dilu_TCarbon(i) * delta_veg(j)) / veget_max(i,j)
281	turnover_daily(i,j,:,:)=(turnover_daily(i,j,:,:)*veget_max_old(i,j)+&
282	dilu_turnover_daily(i,:,:)*delta_veg(j))/veget_max(i,j)
283	bm_to_litter(i,j,:,:)=(bm_to_litter(i,j,:,:)*veget_max_old(i,j)+&
284	dilu_bm_to_litter(i,:,:)*delta_veg(j))/veget_max(i,j)
285	co2flux_new(i,j)=(co2flux_old(i,j)veget_max_old(i,j)+dilu_co2flux_new(i)delta_veg(j))/veget_max(i,j)
286	gpp_daily(i,j)=(gpp_daily(i,j)veget_max_old(i,j)+dilu_gpp_daily(i)delta_veg(j))/veget_max(i,j)
287	resp_growth(i,j)=(resp_growth(i,j)veget_max_old(i,j)+dilu_resp_growth(i)delta_veg(j))/veget_max(i,j)
288	resp_maint(i,j)=(resp_maint(i,j)veget_max_old(i,j)+dilu_resp_maint(i)delta_veg(j))/veget_max(i,j)
289	resp_hetero(i,j)=(resp_hetero(i,j)veget_max_old(i,j)+dilu_resp_hetero(i)delta_veg(j))/veget_max(i,j)
290	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)
291	co2_fire(i,j)=(co2_fire(i,j)veget_max_old(i,j)+dilu_co2_fire(i)delta_veg(j))/veget_max(i,j)
292	ENDIF
293
294	IF(veget_max(i,j).GT.min_stomate) THEN
295
296	! Correct biomass densities to conserve mass
297	! since it's defined on veget_max
298	biomass(i,j,:,:) = biomass(i,j,:,:) * veget_max_old(i,j) / veget_max(i,j)
299
300	ENDIF
301
302	ENDDO ! loop over PFTs
303	ENDDO ! loop over grid points
304
305	vartmp(:)=SUM(gpp_daily*veget_max,dim=2)
306	CALL histwrite_p (hist_id_stomate, "tGPP", itime, vartmp, npts, hori_index)
307	vartmp(:)=SUM(resp_growth*veget_max,dim=2)
308	CALL histwrite_p (hist_id_stomate, "tRESP_GROWTH", itime, vartmp, npts, hori_index)
309	vartmp(:)=SUM(resp_maint*veget_max,dim=2)
310	CALL histwrite_p (hist_id_stomate, "tRESP_MAINT", itime, vartmp, npts, hori_index)
311	vartmp(:)=SUM(resp_hetero*veget_max,dim=2)
312	CALL histwrite_p (hist_id_stomate, "tRESP_HETERO", itime, vartmp, npts, hori_index)
313	vartmp(:)=SUM(co2_to_bm*veget_max,dim=2)
314	CALL histwrite_p (hist_id_stomate, "tCO2_TAKEN", itime, vartmp, npts, hori_index)
315	vartmp(:)=SUM(co2_fire*veget_max,dim=2)
316	CALL histwrite_p (hist_id_stomate, "tCO2_FIRE", itime, vartmp, npts, hori_index)
317	vartmp(:)=SUM(co2flux_new*veget_max,dim=2)
318	CALL histwrite_p (hist_id_stomate, "tCO2FLUX", itime, vartmp, npts, hori_index)
319	vartmp(:)=SUM(co2flux_old*veget_max_old,dim=2)
320	CALL histwrite_p (hist_id_stomate, "tCO2FLUX_OLD", itime, vartmp, npts, hori_index)
321	vartmp(:)=SUM(TCarbon*veget_max,dim=2)
322	CALL histwrite_p (hist_id_stomate, "tCARBON", itime, vartmp, npts, hori_index)
323	vartmp(:)=SUM(SUM(biomass(:,:,:,icarbon),dim=3)*veget_max,dim=2)
324	CALL histwrite_p (hist_id_stomate, "tBIOMASS", itime, vartmp, npts, hori_index)
325	! vartmp(:)=SUM(SUM(SUM(litter(:,:,:,:,icarbon),dim=4),dim=2)*veget_max,dim=2)
326	! CALL histwrite_p (hist_id_stomate, "tLITTER", itime, vartmp, npts, hori_index)
327	! vartmp(:)=SUM(SUM(carbon,dim=2)*veget_max,dim=2)
328	! CALL histwrite_p (hist_id_stomate, "tSOILC", itime, vartmp, npts, hori_index)
329	ENDIF
330
331	END SUBROUTINE cover
332
333	END MODULE lpj_cover

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: branches/publications/ORCHIDEE-LEAK-r5919/src_stomate/lpj_cover.f90 @ 8624

Download in other formats: