source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/stomate_prescribe.f90 @ 6737

Last change on this file since 6737 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 13.3 KB
Line 
1! =================================================================================================================================
2! MODULE       : stomate_prescribe
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         Initialize and update density, crown area.
10!!
11!!\n DESCRIPTION: None
12!!
13!! RECENT CHANGE(S): None
14!!
15!! REFERENCE(S) :
16!!
17!! SVN          :
18!! $HeadURL$
19!! $Date$
20!! $Revision$
21!! \n
22!_ ================================================================================================================================
23
24MODULE stomate_prescribe
25
26  ! modules used:
27
28  USE ioipsl_para
29  USE stomate_data
30  USE pft_parameters
31  USE constantes
32
33  IMPLICIT NONE
34
35  ! private & public routines
36
37  PRIVATE
38  PUBLIC prescribe,prescribe_clear
39
40    ! first call
41    LOGICAL, SAVE                                              :: firstcall_prescribe = .TRUE.
42!$OMP THREADPRIVATE(firstcall_prescribe)
43
44CONTAINS
45
46! =================================================================================================================================
47!! SUBROUTINE   : prescribe_clear
48!!
49!>\BRIEF        : Set the firstcall_prescribe flag back to .TRUE. to prepare for the next simulation.
50!_=================================================================================================================================
51
52  SUBROUTINE prescribe_clear
53    firstcall_prescribe=.TRUE.
54  END SUBROUTINE prescribe_clear
55
56!! ================================================================================================================================
57!! SUBROUTINE   : prescribe
58!!
59!>\BRIEF         Works only with static vegetation and agricultural PFT. Initialize biomass,
60!!               density, crown area in the first call and update them in the following.
61!!
62!! DESCRIPTION (functional, design, flags): \n
63!! This module works only with static vegetation and agricultural PFT.
64!! In the first call, initialize density of individuals, biomass, crown area,
65!! and leaf age distribution to some reasonable value. In the following calls,
66!! these variables are updated.
67!!
68!! To fulfill these purposes, pipe model are used:
69!! \latexonly
70!!     \input{prescribe1.tex}
71!!     \input{prescribe2.tex}
72!! \endlatexonly
73!!
74!! RECENT CHANGE(S): None
75!!
76!! MAIN OUTPUT VARIABLES(S): ::ind, ::cn_ind, ::leaf_frac
77!!
78!! REFERENCES   :
79!! - Krinner, G., N. Viovy, et al. (2005). "A dynamic global vegetation model
80!!   for studies of the coupled atmosphere-biosphere system." Global
81!!   Biogeochemical Cycles 19: GB1015, doi:1010.1029/2003GB002199.
82!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
83!!   plant geography and terrestrial carbon cycling in the LPJ dynamic
84!!   global vegetation model, Global Change Biology, 9, 161-185.
85!!
86!! FLOWCHART    : None
87!! \n
88!_ ================================================================================================================================
89
90 SUBROUTINE prescribe (npts, &
91                        veget_cov_max, dt, PFTpresent, everywhere, when_growthinit, &
92                        biomass, leaf_frac, ind, cn_ind, co2_to_bm)
93
94!! 0. Parameters and variables declaration
95
96   !! 0.1 Input variables
97
98    INTEGER(i_std), INTENT(in)                                :: npts            !! Domain size (unitless)
99    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: veget_cov_max   !! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground (unitless;0-1)
100    REAL(r_std), INTENT(in)                                   :: dt              !! time step (dt_days)
101    !! 0.2 Output variables
102
103    !! 0.3 Modified variables
104
105    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)               :: PFTpresent      !! PFT present (0 or 1)
106    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: everywhere      !! is the PFT everywhere in the grid box or very localized (after its introduction) (?)
107    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: when_growthinit !! how many days ago was the beginning of the growing season (days)
108    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass   !! biomass (gC/(m^2 of ground))
109    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac       !! fraction of leaves in leaf age class (unitless;0-1)
110    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: ind             !! density of individuals (1/(m^2 of ground))
111    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: cn_ind          !! crown area per individual (m^2)
112    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: co2_to_bm       !! co2 taken up by carbohydrate
113                                                                                 !! reserve at the beginning of the
114                                                                                 !! growing season @tex ($gC m^{-2}
115                                                                                 !! of total ground/day$) @endtex
116
117    !! 0.4 Local variables
118
119    REAL(r_std), DIMENSION(npts)                              :: dia             !! stem diameter (m)
120    REAL(r_std), DIMENSION(npts)                              :: woodmass        !! woodmass (gC/(m^2 of ground))
121    REAL(r_std), DIMENSION(npts)                              :: woodmass_ind    !! woodmass of an individual (gC)
122    INTEGER(i_std)                                            :: i,j             !! index (unitless)
123
124!_ ================================================================================================================================
125
126    DO j = 2,nvm
127
128      ! only when the DGVM is not activated or agricultural PFT.
129
130      IF ( ( .NOT. ok_dgvm .AND. lpj_gap_const_mort ) .OR. ( .NOT. natural(j) .OR. pasture(j)) ) THEN
131
132        !
133        !! 1.Update crown area
134        !
135
136        cn_ind(:,j) = zero
137
138        IF ( is_tree(j) ) THEN
139
140          !
141          !! 1.1 treat for trees
142          !
143
144          dia(:) = zero
145
146          DO i = 1, npts ! loop over grid points
147
148            IF ( veget_cov_max(i,j) .GT. zero ) THEN
149
150              !! 1.1.1 calculate wood mass on an area basis, which include sapwood and heartwood aboveground and belowground.
151
152              woodmass(i) = (biomass(i,j,isapabove,icarbon) + biomass(i,j,isapbelow,icarbon) + &
153                   biomass(i,j,iheartabove,icarbon) + biomass(i,j,iheartbelow,icarbon)) * veget_cov_max(i,j) 
154
155              IF ( woodmass(i) .GT. min_stomate ) THEN
156
157                !! 1.1.2 calculate critical individual density
158!?? the logic for 1.1.3 and 1.1.2 is strange, it should be the case that first to calculate critical woodmass per individual,
159!?? then calculate critical density.
160
161
162                ! how to derive the following equation:
163                ! first, TreeHeight=pipe_tune2 * Diameter^{pipe_tune3}
164                ! we assume the tree is an ideal cylinder, so it volume is: Volume = pi*(Dia/2)^2*Height = pi/4 * Dia * pipe_tune2*Dia^{pipe_tune3}
165                !                                                                  = pi/4 * pipe_tune2 * Dia^{2+pipe_tune3}
166                ! last, the woodmass_per_individual = pipe_density * Volume = pipe_density*pi/4.*pipe_tune2 * Dia^{2+pipe_tune3}             
167                ind(i,j) = woodmass(i) / &
168                           ( pipe_density*pi/4.*pipe_tune2 * maxdia(j)**(2.+pipe_tune3) )
169
170                !! 1.1.3 individual biomass corresponding to this critical density of individuals
171
172                woodmass_ind(i) = woodmass(i) / ind(i,j)
173
174                !! 1.1.4 calculate stem diameter per individual tree
175
176                dia(i) = ( woodmass_ind(i) / ( pipe_density * pi/4. * pipe_tune2 ) ) ** &
177                         ( un / ( 2. + pipe_tune3 ) )
178
179                !! 1.1.5 calculate provisional tree crown area for per individual tree
180
181                ! equation: CrownArea=pipe_tune1 * Diameter^{1.6}
182                cn_ind(i,j) = pipe_tune1 * MIN( maxdia(j), dia(i) ) ** pipe_tune_exp_coeff
183
184                !! 1.1.6 If total tree crown area for this tree PFT exceeds its veget_cov_max, tree density is recalculated.
185
186                IF ( cn_ind(i,j) * ind(i,j) .GT. 1.002* veget_cov_max(i,j) ) THEN
187
188                  ind(i,j) = veget_cov_max(i,j) / cn_ind(i,j)
189
190                ELSE
191
192                   ind(i,j) = ( veget_cov_max(i,j) / &
193                        &     ( pipe_tune1 * (woodmass(i)/(pipe_density*pi/4.*pipe_tune2)) &
194                        &     **(pipe_tune_exp_coeff/(2.+pipe_tune3)) ) ) &
195                        &     ** (1./(1.-(pipe_tune_exp_coeff/(2.+pipe_tune3))))
196                   
197
198                  woodmass_ind(i) = woodmass(i) / ind(i,j)
199
200                  dia(i) = ( woodmass_ind(i) / ( pipe_density * pi/4. * pipe_tune2 ) ) ** &
201                           ( un / ( 2. + pipe_tune3 ) )
202
203                  ! final crown area
204                  cn_ind(i,j) = pipe_tune1 * MIN( maxdia(j), dia(i) ) ** pipe_tune_exp_coeff
205
206                ENDIF
207
208              ELSE !woodmas=0  => impose some value
209
210                dia(:) = maxdia(j)
211
212                cn_ind(i,j) = pipe_tune1 * MIN( maxdia(j), dia(i) ) ** pipe_tune_exp_coeff
213
214              ENDIF ! IF ( woodmass(i) .GT. min_stomate )
215
216            ENDIF    ! veget_cov_max .GT. 0.
217
218          ENDDO      ! loop over grid points
219
220        ELSE !grass
221
222          !
223          !! 1.2 grasses: crown area always set to 1m**2
224          !
225
226          WHERE ( veget_cov_max(:,j) .GT. zero )
227            cn_ind(:,j) = un
228          ENDWHERE
229
230        ENDIF   !IF ( is_tree(j) )
231
232        !
233        !! 2 density of individuals
234        !
235       
236        WHERE ( veget_cov_max(:,j) .GT. zero )
237
238          ind(:,j) = veget_cov_max(:,j) / cn_ind(:,j) 
239
240        ELSEWHERE
241
242          ind(:,j) = zero
243
244        ENDWHERE
245
246      ENDIF     ! IF ( ( .NOT. ok_dgvm .AND. lpj_gap_const_mort ) .OR. ( .NOT. natural(j) ) )
247
248    ENDDO       ! loop over PFTs
249
250    !
251    !!? it's better to move the code for first call at the beginning of the module.
252    !! 2 If it's the first call for this module,
253    !
254
255    IF (( firstcall_prescribe ) .AND. (TRIM(stom_restname_in) == 'NONE')) THEN
256
257       IF (printlev >= 2) THEN
258          WRITE(numout,*) 'prescribe:'
259          ! impose some biomass if zero and PFT prescribed
260          WRITE(numout,*) '   > Imposing initial biomass for prescribed trees, '// &
261               'initial reserve mass for prescribed grasses.'
262          WRITE(numout,*) '   > Declaring prescribed PFTs present.'
263       END IF
264
265      DO j = 2,nvm ! loop over PFTs
266        DO i = 1, npts ! loop over grid points
267
268          ! is vegetation static or PFT agricultural?
269          ! Static vegetation or agricultural PFT
270          IF ( ( .NOT. ok_dgvm ) .OR. &
271               ( ( .NOT. natural(j) .OR. pasture(j)) .AND. ( veget_cov_max(i,j) .GT. min_stomate ) ) ) THEN
272
273            !
274            !! 2.1 if tree biomass is extremely small, prescribe the biomass by assuming they have sapling biomass, which is a constant in the model.
275            !!     then set all the leaf age as 1.
276            !
277            ! if tree PFT and biomass too small, prescribe the biomass to a value.
278            IF ( is_tree(j) .AND. &
279                 ( veget_cov_max(i,j) .GT. min_stomate ) .AND. &
280                 ( SUM( biomass(i,j,:,icarbon) ) .LE. min_stomate ) ) THEN
281               !!? here the code is redundant, as "veget_cov_max(i,j) .GT. min_stomate" is already met in the above if condition.
282               IF (veget_cov_max(i,j) .GT. min_stomate) THEN
283                  biomass(i,j,:,:) = (bm_sapl_rescale * bm_sapl(j,:,:) * ind(i,j)) / veget_cov_max(i,j)
284               ELSE
285                  biomass(i,j,:,:) = zero
286               ENDIF
287
288              ! set leaf age classes
289              leaf_frac(i,j,:) = zero
290              leaf_frac(i,j,1) = un
291
292              ! set time since last beginning of growing season
293              when_growthinit(i,j) = large_value
294
295              ! seasonal trees: no leaves at beginning
296
297              IF ( pheno_model(j) .NE. 'none' ) THEN
298
299                biomass(i,j,ileaf,icarbon) = zero
300                leaf_frac(i,j,1) = zero
301
302              ENDIF
303
304              co2_to_bm(i,j) = co2_to_bm(i,j) + ( SUM(biomass(i,j,:,icarbon))  / dt )
305            ENDIF
306
307            !
308            !! 2.2 for grasses, set only the carbon reserve pool to "sapling" carbon reserve pool.
309            !!     and set all leaf age to 1.
310
311            IF ( ( .NOT. is_tree(j) ) .AND. &
312                 ( veget_cov_max(i,j) .GT. min_stomate ) .AND. &
313                 ( SUM( biomass(i,j,:,icarbon) ) .LE. min_stomate ) ) THEN
314
315              biomass(i,j,icarbres,:) = bm_sapl(j,icarbres,:) * ind(i,j) / veget_cov_max(i,j)
316
317              ! set leaf age classes
318              leaf_frac(i,j,:) = zero
319              leaf_frac(i,j,1) = un
320
321              ! set time since last beginning of growing season
322              when_growthinit(i,j) = large_value
323
324              co2_to_bm(i,j) = co2_to_bm(i,j) + ( biomass(i,j,icarbres,icarbon)  / dt )
325            ENDIF
326
327            !
328            !! 2.3 declare all PFTs with positive veget_cov_max as present everywhere in that grid box
329            !
330
331            IF ( veget_cov_max(i,j) .GT. min_stomate ) THEN
332              PFTpresent(i,j) = .TRUE.
333              everywhere(i,j) = un
334            ENDIF
335
336          ENDIF   ! not ok_dgvm  or agricultural
337
338        ENDDO ! loop over grid points
339      ENDDO ! loop over PFTs
340
341
342   ENDIF
343   
344   firstcall_prescribe = .FALSE.
345
346  END SUBROUTINE prescribe
347
348END MODULE stomate_prescribe
Note: See TracBrowser for help on using the repository browser.