| 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, height, basal area and 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 | |
|---|
| 24 | MODULE stomate_prescribe |
|---|
| 25 | |
|---|
| 26 | ! modules used: |
|---|
| 27 | |
|---|
| 28 | USE ioipsl_para |
|---|
| 29 | USE xios_orchidee |
|---|
| 30 | USE time, ONLY : LastTsYear |
|---|
| 31 | USE pft_parameters |
|---|
| 32 | USE constantes |
|---|
| 33 | USE function_library, ONLY: calculate_c0_alloc, cc_to_lai, & |
|---|
| 34 | cc_to_biomass, check_vegetation_area, check_mass_balance, & |
|---|
| 35 | wood_to_qmdia, Nmax, sort_circ_class_biomass, lai_to_biomass, & |
|---|
| 36 | biomass_to_lai, wood_to_ba |
|---|
| 37 | |
|---|
| 38 | IMPLICIT NONE |
|---|
| 39 | |
|---|
| 40 | ! private & public routines |
|---|
| 41 | |
|---|
| 42 | PRIVATE |
|---|
| 43 | PUBLIC prescribe, prescribe_clear |
|---|
| 44 | |
|---|
| 45 | ! first call |
|---|
| 46 | LOGICAL, SAVE :: firstcall_prescribe = .TRUE. |
|---|
| 47 | !$OMP THREADPRIVATE(firstcall_prescribe) |
|---|
| 48 | INTEGER(i_std), SAVE :: printlev_loc !! Local level of text output for this module |
|---|
| 49 | !$OMP THREADPRIVATE(printlev_loc) |
|---|
| 50 | |
|---|
| 51 | CONTAINS |
|---|
| 52 | |
|---|
| 53 | ! ================================================================================================================================= |
|---|
| 54 | !! SUBROUTINE : prescribe_clear |
|---|
| 55 | !! |
|---|
| 56 | !>\BRIEF : Set the firstcall_prescribe flag back to .TRUE. to prepare for the next simulation. |
|---|
| 57 | !_================================================================================================================================= |
|---|
| 58 | |
|---|
| 59 | SUBROUTINE prescribe_clear |
|---|
| 60 | firstcall_prescribe=.TRUE. |
|---|
| 61 | END SUBROUTINE prescribe_clear |
|---|
| 62 | |
|---|
| 63 | |
|---|
| 64 | !! ================================================================================================================================ |
|---|
| 65 | !! SUBROUTINE : prescribe |
|---|
| 66 | !! |
|---|
| 67 | !>\BRIEF Works only with static vegetation and agricultural PFTs. |
|---|
| 68 | !! Initialize biomass, density, presence in the first call |
|---|
| 69 | !! and update them in the following. |
|---|
| 70 | !! |
|---|
| 71 | !! DESCRIPTION (functional, design, flags): \n |
|---|
| 72 | !! This module works only with static vegetation and agricultural PFTs. |
|---|
| 73 | !! In the first call, initialize density of individuals, biomass, |
|---|
| 74 | !! and leaf age distribution to some reasonable value. In the following calls, |
|---|
| 75 | !! these variables are updated following a stand replacing disturbance or |
|---|
| 76 | !! if the conditions are suitable for recruitment. |
|---|
| 77 | !! |
|---|
| 78 | !! To fulfill these purposes, pipe model are used: |
|---|
| 79 | !! \latexonly |
|---|
| 80 | !! \input{prescribe1.tex} |
|---|
| 81 | !! \input{prescribe2.tex} |
|---|
| 82 | !! \endlatexonly |
|---|
| 83 | !! |
|---|
| 84 | !! RECENT CHANGE(S): None |
|---|
| 85 | !! |
|---|
| 86 | !! MAIN OUTPUT VARIABLES(S): ::ind, ::biomass |
|---|
| 87 | !! |
|---|
| 88 | !! REFERENCES : |
|---|
| 89 | !! - Krinner, G., N. Viovy, et al. (2005). "A dynamic global vegetation model |
|---|
| 90 | !! for studies of the coupled atmosphere-biosphere system." Global |
|---|
| 91 | !! Biogeochemical Cycles 19: GB1015, doi:1010.1029/2003GB002199. |
|---|
| 92 | !! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics, |
|---|
| 93 | !! plant geography and terrestrial carbon cycling in the LPJ dynamic |
|---|
| 94 | !! global vegetation model, Global Change Biology, 9, 161-185. |
|---|
| 95 | !! - McDowell, N., Barnard, H., Bond, B.J., Hinckley, T., Hubbard, R.M., Ishii, |
|---|
| 96 | !! H., Köstner, B., Magnani, F. Marshall, J.D., Meinzer, F.C., Phillips, N., |
|---|
| 97 | !! Ryan, M.G., Whitehead D. 2002. The relationship between tree height and leaf |
|---|
| 98 | !! area: sapwood area ratio. Oecologia, 132:12â20. |
|---|
| 99 | !! - Novick, K., Oren, R., Stoy, P., Juang, F.-Y., Siqueira, M., Katul, G. 2009. |
|---|
| 100 | !! The relationship between reference canopy conductance and simplified hydraulic |
|---|
| 101 | !! architecture. Advances in water resources 32, 809-819. |
|---|
| 102 | !! - Ruger et al. 2009, J. of Ecol. doi: 10.1111/j.1365-2745.2009.01552.x |
|---|
| 103 | !! |
|---|
| 104 | !! FLOWCHART : None |
|---|
| 105 | !! \n |
|---|
| 106 | !_ ================================================================================================================================ |
|---|
| 107 | |
|---|
| 108 | SUBROUTINE prescribe (& |
|---|
| 109 | npts, veget_max, dt, PFTpresent, & |
|---|
| 110 | everywhere, when_growthinit, leaf_frac, circ_class_dist, & |
|---|
| 111 | circ_class_n, circ_class_biomass, atm_to_bm, forest_managed, & |
|---|
| 112 | KF, plant_status, age, npp_longterm, & |
|---|
| 113 | lm_lastyearmax, longevity_eff_leaf, longevity_eff_sap, longevity_eff_root, & |
|---|
| 114 | k_latosa_adapt, light_tran_to_floor_season, lpft_replant, species_change_map, & |
|---|
| 115 | cn_leaf_init_2D, bm_sapl_2D, new_ind, qmd_init, & |
|---|
| 116 | dia_init) |
|---|
| 117 | |
|---|
| 118 | !! 0. Parameters and variables declaration |
|---|
| 119 | |
|---|
| 120 | !! 0.1 Input variables |
|---|
| 121 | |
|---|
| 122 | INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless) |
|---|
| 123 | INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: forest_managed !! forest management flag |
|---|
| 124 | REAL(r_std), INTENT(in) :: dt !! time step (dt_days) |
|---|
| 125 | REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT |
|---|
| 126 | !! (LAI -> infinity) on ground. May sum to |
|---|
| 127 | !! less than unity if the pixel has |
|---|
| 128 | !! nobio area. (unitless; 0-1) |
|---|
| 129 | REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: longevity_eff_root !! Effective root turnover time that accounts |
|---|
| 130 | !! waterstress (days) |
|---|
| 131 | REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: longevity_eff_sap !! Effective sapwood turnover time that accounts |
|---|
| 132 | !! waterstress (days) |
|---|
| 133 | REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: longevity_eff_leaf !! Effective leaf turnover time that accounts |
|---|
| 134 | !! waterstress (days) |
|---|
| 135 | LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: lpft_replant !! Set to true if a PFT has been clearcut |
|---|
| 136 | !! and needs to be replaced by another species |
|---|
| 137 | INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: species_change_map !! map which gives the PFT number that each PFT |
|---|
| 138 | REAL(r_std), DIMENSION(npts,nvm),INTENT(inout) :: light_tran_to_floor_season !! Seasonal fraction of light transmitted to canopy levels |
|---|
| 139 | REAL(r_std), DIMENSION(ncirc), INTENT(in) :: circ_class_dist !! The probability distribution of trees |
|---|
| 140 | !! in a circ class in case of a |
|---|
| 141 | !! redistribution (unitless). |
|---|
| 142 | REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: cn_leaf_init_2D !! initial leaf C/N ratio |
|---|
| 143 | REAL(r_std), DIMENSION(nvm), INTENT(in) :: qmd_init !! quadratic mean diameter of a newly planted PFT (m) |
|---|
| 144 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: dia_init !! Initial diameter distribution of a newly planted PFT (m) |
|---|
| 145 | |
|---|
| 146 | !! 0.2 Output variables |
|---|
| 147 | REAL(r_std),DIMENSION(npts,nvm), INTENT(out) :: new_ind !! density of new individuals (trees m-2 day-1) |
|---|
| 148 | |
|---|
| 149 | !! 0.3 Modified variables |
|---|
| 150 | |
|---|
| 151 | LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: PFTpresent !! PFT present (0 or 1) |
|---|
| 152 | REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: plant_status !! Growth and phenological status of the plant |
|---|
| 153 | !! Different stati are listed in in constantes_var |
|---|
| 154 | REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or |
|---|
| 155 | !! very localized (after its introduction) (?) |
|---|
| 156 | REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: when_growthinit !! how many days ago was the beginning of |
|---|
| 157 | !! the growing season (days) |
|---|
| 158 | REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: npp_longterm !! "long term" net primary productivity |
|---|
| 159 | !! @tex ($gC m^{-2} year^{-1}$) @endtex |
|---|
| 160 | REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age !! mean age (years) |
|---|
| 161 | REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT |
|---|
| 162 | !! @tex ($gC m^{-2}$) @endtex |
|---|
| 163 | REAL(r_std), DIMENSION(npts,nvm,ncirc), INTENT(inout) :: circ_class_n !! Number of individuals in each circ class |
|---|
| 164 | !! @tex $(ind m^{-2})$ @endtex |
|---|
| 165 | REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), INTENT(inout) :: circ_class_biomass !! Biomass components of the model tree |
|---|
| 166 | !! within a circumference class |
|---|
| 167 | !! class @tex $(g C ind^{-1})$ @endtex |
|---|
| 168 | REAL(r_std), DIMENSION(npts,nvm,nelements), INTENT(inout) :: atm_to_bm !! CO2 or N taken from the atmosphere to get |
|---|
| 169 | !! the seed C and N |
|---|
| 170 | !! to create the seedlings |
|---|
| 171 | !! @tex (gC.m^{-2}dt^{-1})$ @endtex |
|---|
| 172 | REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age |
|---|
| 173 | !! class (unitless;0-1) |
|---|
| 174 | REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: KF !! Scaling factor to convert sapwood mass |
|---|
| 175 | !! into leaf mass (m) - this variable is |
|---|
| 176 | !! passed to other routines |
|---|
| 177 | REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: k_latosa_adapt !! Leaf to sapwood area adapted for long |
|---|
| 178 | !! term water stress (m) |
|---|
| 179 | REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), INTENT(inout) :: bm_sapl_2D !! Sapling biomass for the functional |
|---|
| 180 | !! allocation with a dimension for the |
|---|
| 181 | !! circumference classes |
|---|
| 182 | !! @tex $(gC.ind^{-1})$ @endtex |
|---|
| 183 | |
|---|
| 184 | !! 0.4 Local variables |
|---|
| 185 | REAL(r_std), DIMENSION(nvm) :: c0_alloc !! Root to sapwood tradeoff parameter |
|---|
| 186 | INTEGER(i_std) :: ipts, ivm !! index (unitless) |
|---|
| 187 | INTEGER(i_std) :: ipar, ilev !! index (unitless) |
|---|
| 188 | INTEGER(i_std) :: icir, iele, imbc !! index (unitless) |
|---|
| 189 | INTEGER(i_std) :: deb,fin, imaxt !! index (unitless) |
|---|
| 190 | REAL(r_std), DIMENSION(npts,nvm) :: LF !! Scaling factor to convert sapwood mass |
|---|
| 191 | !! into root mass (unitless) |
|---|
| 192 | REAL(r_std), DIMENSION(npts,nvm) :: lstress_fac !! Light stress factor, based on total |
|---|
| 193 | !! transmitted light (unitless, 0-1) |
|---|
| 194 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: circ_class !! circumference of individual trees |
|---|
| 195 | REAL(r_std), DIMENSION(nparts) :: bm_init !! Biomass needed to initiate the next |
|---|
| 196 | !! planting @tex $(gC m^{-2})$ @endtex |
|---|
| 197 | REAL(r_std), DIMENSION(ncirc) :: nb_trees_i !! Number of trees in each twentith |
|---|
| 198 | !! circumference quantile of the |
|---|
| 199 | !! distribution (ind) |
|---|
| 200 | REAL(r_std) :: excedent !! Number of trees after truncation to be |
|---|
| 201 | !! reallocated to smaller quantiles of the |
|---|
| 202 | !! distribution (ind) |
|---|
| 203 | REAL(r_std) :: max_circ_init !! Maximum initial circumferences of the |
|---|
| 204 | !! truncated exponential distribution (cm) |
|---|
| 205 | REAL(r_std), DIMENSION(ncirc) :: ave_tree_height !! The height of the ideal tree in each |
|---|
| 206 | REAL(r_std) :: dia_init_recr !! Initial diameter of a recruit (m). Note that |
|---|
| 207 | !! this value may differ from the input argument |
|---|
| 208 | !! dia_int which is the diameter distribution of |
|---|
| 209 | !! a newly planted PFT |
|---|
| 210 | REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements):: check_intern !! Contains the components of the internal |
|---|
| 211 | !! mass balance chech for this routine |
|---|
| 212 | !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex |
|---|
| 213 | REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance |
|---|
| 214 | !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex |
|---|
| 215 | REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start and end pool of this routine |
|---|
| 216 | !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex |
|---|
| 217 | REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! Start and end pool of this routine |
|---|
| 218 | !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex |
|---|
| 219 | REAL(r_std) :: sla_est !! A first estimate of sla in case its calculation is |
|---|
| 220 | !! dynamic @tex $(m^2.gC^{-1})$ @endtex |
|---|
| 221 | REAL(r_std), DIMENSION(ncirc) :: bypass_lm !! leaf mass before it is set to zero |
|---|
| 222 | !! for deciduous PFTs. This value is |
|---|
| 223 | !! in the calculation of KF (gC tree-1) |
|---|
| 224 | REAL(r_std), DIMENSION(npts,nvm) :: veget_max_begin !! temporary storage of veget_max to check area conservation |
|---|
| 225 | REAL(r_std) :: sapwood_density !! Sapwood density gC m^{-3} |
|---|
| 226 | REAL(r_std) :: ind_init !! Maximum number of individuals when establishing a new forest (#/m^{-2}) |
|---|
| 227 | REAL(r_std) :: wood_init !! Wood mass of a sapling (g C tree^{-1})i |
|---|
| 228 | REAL(r_std),DIMENSION(ncirc) :: ba_init !! intitial basal area |
|---|
| 229 | LOGICAL :: establish !! Logical flag to establish a new young vegetation |
|---|
| 230 | !! after a stand replacing disturbance |
|---|
| 231 | LOGICAL :: recruite !! Logical flag to add saplings of a PFT under an |
|---|
| 232 | !! existing canopy of the same PFT |
|---|
| 233 | REAL(r_std), DIMENSION(npts,nvm,ncirc) :: old_circ_class_n !! Old number of individuals in each circ class |
|---|
| 234 | !! @tex $(ind m^{-2})$ @endtex |
|---|
| 235 | REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements) & |
|---|
| 236 | :: old_circ_class_biomass !! Old biomass components of the model tree |
|---|
| 237 | !! within a circumference class |
|---|
| 238 | !! class @tex $(g C ind^{-1})$ @endtex |
|---|
| 239 | REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: old_biomass !! Old stand level biomass |
|---|
| 240 | !! tex $(gC.m^{-2})$ @endtex |
|---|
| 241 | REAL(r_std), DIMENSION(ncirc) :: height_small !! Height of the smallest tree in circ_class_biomass |
|---|
| 242 | !! tex $(m)$ @endtex |
|---|
| 243 | REAL(r_std) :: dia_small !! Diameter of the smallest tree in circ_class_biomass |
|---|
| 244 | !! tex $(m)$ @endtex |
|---|
| 245 | REAL(r_std) :: height_sapl !! Height of the recruitment sapling |
|---|
| 246 | !! tex $(m)$ @endtex |
|---|
| 247 | REAL(r_std) :: old_ind !! density of old individuals |
|---|
| 248 | REAL(r_std) :: cn_leaf,cn_root !! CN ratio of leaves, and roots |
|---|
| 249 | !! (gC/gN) |
|---|
| 250 | REAL(r_std) :: cn_wood !! CN ratio of wood pool |
|---|
| 251 | !! (gC/gN) |
|---|
| 252 | REAL(r_std) :: coeff |
|---|
| 253 | REAL(r_std), DIMENSION(ncirc) :: dummy |
|---|
| 254 | |
|---|
| 255 | REAL(r_std) :: lambda !! lambda of the truncated exponential |
|---|
| 256 | !! distribution of initial diameters (1/cm**2) |
|---|
| 257 | REAL(r_std) :: tmp_ind !! temporary number of individuals per m-2 |
|---|
| 258 | REAL(r_std), DIMENSION(ncirc) :: circ_class_n_new !! variable used to sort circ_class_n |
|---|
| 259 | REAL(r_std), DIMENSION(ncirc,nparts,nelements) :: circ_class_biomass_new !! variable used to sort circ_class_biomass |
|---|
| 260 | REAL(r_std) :: k_latosa_tmp !! Temporaray variable in the calculation of |
|---|
| 261 | !! k_latosa_adapt |
|---|
| 262 | REAL(r_std), DIMENSION(nvm) :: n_init !! Initial number of individuals calculated |
|---|
| 263 | !! based on gradient calculation (ind m^(-2)) |
|---|
| 264 | REAL(r_std), DIMENSION(npts,nvm) :: tmp !! Temporary variable to prepare information for xios |
|---|
| 265 | !_ ================================================================================================================================ |
|---|
| 266 | |
|---|
| 267 | |
|---|
| 268 | IF (firstcall_prescribe) THEN |
|---|
| 269 | !! Initialize local printlev |
|---|
| 270 | printlev_loc=get_printlev('stomate_prescribe') |
|---|
| 271 | firstcall_prescribe=.FALSE. |
|---|
| 272 | END IF |
|---|
| 273 | |
|---|
| 274 | IF (printlev_loc>=2) WRITE(numout,*) 'Entering stomate_prescribe' |
|---|
| 275 | |
|---|
| 276 | !! 1. Initialize biomass at first call |
|---|
| 277 | |
|---|
| 278 | !! 1.0 Initialize variables that won't be calculated for all PFTs, |
|---|
| 279 | !! to avoid a NaN if the value for bare soil is ever written. |
|---|
| 280 | c0_alloc(:)=zero |
|---|
| 281 | |
|---|
| 282 | !! 1.1 Store for closing the mass balance with recruitment |
|---|
| 283 | old_biomass(:,:,:,:)= cc_to_biomass(npts,nvm,& |
|---|
| 284 | circ_class_biomass(:,:,:,:,:), & |
|---|
| 285 | circ_class_n(:,:,:)) |
|---|
| 286 | new_ind(:,:) = zero |
|---|
| 287 | |
|---|
| 288 | !! 1.3 Initialize check for mass balance closure |
|---|
| 289 | IF (err_act.GT.1) THEN |
|---|
| 290 | |
|---|
| 291 | check_intern(:,:,:,:) = zero |
|---|
| 292 | pool_start(:,:,:) = zero |
|---|
| 293 | DO iele = 1,nelements |
|---|
| 294 | |
|---|
| 295 | ! atm_to_bm has as intent inout, the variable |
|---|
| 296 | ! accumulates carbon over the course of a day. |
|---|
| 297 | ! Use the difference between start and the end of |
|---|
| 298 | ! this routine |
|---|
| 299 | check_intern(:,:,iatm2land,iele) = - un * & |
|---|
| 300 | atm_to_bm(:,:,iele) * veget_max(:,:) * dt |
|---|
| 301 | |
|---|
| 302 | DO ipar = 1,nparts |
|---|
| 303 | |
|---|
| 304 | DO icir = 1,ncirc |
|---|
| 305 | |
|---|
| 306 | ! Initial biomass pool. Use circ_class_biomass |
|---|
| 307 | ! because that variable is the prognostic variable |
|---|
| 308 | ! biomass is syncronized to circ_class_biomass. If |
|---|
| 309 | ! many diameter classes are used, each diameter class |
|---|
| 310 | ! separatly passes all criteria but when combined |
|---|
| 311 | ! a mass balance problem occurs. |
|---|
| 312 | pool_start(:,:,iele) = pool_start(:,:,iele) + & |
|---|
| 313 | circ_class_biomass(:,:,icir,ipar,iele) * & |
|---|
| 314 | circ_class_n(:,:,icir) * veget_max(:,:) |
|---|
| 315 | |
|---|
| 316 | ENDDO |
|---|
| 317 | |
|---|
| 318 | ENDDO |
|---|
| 319 | |
|---|
| 320 | ENDDO |
|---|
| 321 | |
|---|
| 322 | !! 1.3 Initialize check for area conservation |
|---|
| 323 | veget_max_begin(:,:) = veget_max(:,:) |
|---|
| 324 | |
|---|
| 325 | ENDIF ! err_act.GT.1 |
|---|
| 326 | |
|---|
| 327 | !! 2. Prescribe carbon and nitrogen biomass |
|---|
| 328 | |
|---|
| 329 | ! We want the vegetation to regrow after a stand replacing disturbance |
|---|
| 330 | ! So loop over every point and PFT. If there is supposed to be |
|---|
| 331 | ! vegetation here (according to veget_max) but there isn't (according |
|---|
| 332 | ! to ind), then establish a new young vegetation in that PFT. |
|---|
| 333 | DO ivm = 2,nvm ! Loop over PFTs |
|---|
| 334 | |
|---|
| 335 | DO ipts = 1, npts ! Loop over pixels |
|---|
| 336 | |
|---|
| 337 | !! 2.1 Wait for the end of the year to replant |
|---|
| 338 | ! If we are regrowing different species after managed forests |
|---|
| 339 | ! are killed/die, we sometimes need to prevent them regrowing. |
|---|
| 340 | ! This is only true if they die in the middle of the year due |
|---|
| 341 | ! to natural causes, in which case they will be replanted at |
|---|
| 342 | ! the end of the year. This loop checks to see if we regrow |
|---|
| 343 | ! this PFT now or later. We can regrow a PFT after lpft_replant |
|---|
| 344 | ! is set back to false so before that we need to process all the |
|---|
| 345 | ! information such that the correct species will be regrown. |
|---|
| 346 | ! lpft_replant is an optional variable. First test whether it is |
|---|
| 347 | ! present. |
|---|
| 348 | IF(ok_change_species)THEN |
|---|
| 349 | IF(lpft_replant(ipts,ivm))THEN |
|---|
| 350 | CYCLE |
|---|
| 351 | ENDIF |
|---|
| 352 | ENDIF |
|---|
| 353 | |
|---|
| 354 | ! Debug |
|---|
| 355 | IF(printlev_loc.GE.4 .AND. ipts == test_grid & |
|---|
| 356 | .AND. ivm == test_pft)THEN |
|---|
| 357 | WRITE(numout,*) 'PFT ',ivm |
|---|
| 358 | WRITE(numout,*) 'plant_status ', plant_status(ipts,ivm) |
|---|
| 359 | WRITE(numout,*) 'veget_max ', veget_max(ipts,ivm) |
|---|
| 360 | WRITE(numout,*) 'recruitment_pft ',recruitment_pft(ivm) |
|---|
| 361 | WRITE(numout,*) 'ok_dgvm ', ok_dgvm |
|---|
| 362 | WRITE(numout,*) 'biomass - C, ', & |
|---|
| 363 | SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) |
|---|
| 364 | WRITE(numout,*) 'biomass - N, ', & |
|---|
| 365 | SUM(SUM(circ_class_biomass(ipts,ivm,:,:,initrogen),1)) |
|---|
| 366 | WRITE(numout,*) 'LastTsYear ', LastTsYear |
|---|
| 367 | WRITE(numout,*) 'circ_class_n, ', SUM(circ_class_n(ipts,ivm,:)) |
|---|
| 368 | WRITE(numout,*) 'dens_target, ', dens_target(ivm)*ha_to_m2 |
|---|
| 369 | WRITE(numout,*) 'do_now_recruit, ',do_now_recruit |
|---|
| 370 | ENDIF |
|---|
| 371 | ! |
|---|
| 372 | |
|---|
| 373 | ! Decide whether we have to establish new vegetation (::establish) |
|---|
| 374 | ! allow recruitment (::recruite). These two local variables were |
|---|
| 375 | ! introduced to allow maximal flexibility so within a single |
|---|
| 376 | ! run, PFTs with and without recruitment can co-exist. |
|---|
| 377 | IF( ( (plant_status(ipts,ivm) .EQ. iprescribe) .OR. & |
|---|
| 378 | (plant_status(ipts,ivm) .EQ. idead) ) .AND. & |
|---|
| 379 | (SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .LE. min_stomate) .AND. & |
|---|
| 380 | (veget_max(ipts,ivm) .GT. min_stomate) .AND. & |
|---|
| 381 | (.NOT. ok_dgvm) ) THEN |
|---|
| 382 | |
|---|
| 383 | ! The vegetation just died from replacing disturbances |
|---|
| 384 | ! (two first conditions). The new vegetation could be |
|---|
| 385 | ! of the same or a different PFT than the previous |
|---|
| 386 | ! vegetation. |
|---|
| 387 | establish = .TRUE. |
|---|
| 388 | recruite = .FALSE. |
|---|
| 389 | |
|---|
| 390 | ! Debug |
|---|
| 391 | IF (printlev_loc>=4) THEN |
|---|
| 392 | WRITE(numout,*) 'We are calling stomate_prescribe.f90 & |
|---|
| 393 | & to initialize PFT ',ivm |
|---|
| 394 | ENDIF |
|---|
| 395 | !- |
|---|
| 396 | |
|---|
| 397 | ELSEIF( ( (plant_status(ipts,ivm) .NE. iprescribe) .OR. & |
|---|
| 398 | (plant_status(ipts,ivm) .NE. idead) ) .AND. & |
|---|
| 399 | (SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .GT. min_stomate) .AND. & |
|---|
| 400 | (veget_max(ipts,ivm) .GT. min_stomate) .AND. & |
|---|
| 401 | (recruitment_pft(ivm)) .AND. & |
|---|
| 402 | (.NOT. ok_dgvm) .AND. & |
|---|
| 403 | (LastTsYear) .AND. & |
|---|
| 404 | (SUM(circ_class_n(ipts,ivm,:)) .GT. dens_target(ivm)*ha_to_m2) .AND. & |
|---|
| 405 | do_now_recruit) THEN |
|---|
| 406 | |
|---|
| 407 | ! There is already vegetation (first three conditions) |
|---|
| 408 | ! and this specific PFT has recruitment (condition 5) |
|---|
| 409 | ! Only use recruitment when the recruitment flag is TRUE |
|---|
| 410 | ! (condition 7) |
|---|
| 411 | ! no recruit when the stand density is below dens_target because |
|---|
| 412 | ! the model will slowly but surely kill the forest before a new |
|---|
| 413 | ! establishment |
|---|
| 414 | establish = .FALSE. |
|---|
| 415 | recruite = .TRUE. |
|---|
| 416 | |
|---|
| 417 | ! Debug |
|---|
| 418 | IF (printlev_loc>=4) THEN |
|---|
| 419 | WRITE(numout,*) 'We are calling stomate_prescribe.f90 & |
|---|
| 420 | & for sapling recruitment in PFT ',ivm |
|---|
| 421 | WRITE(numout,*) 'Density of stem', SUM(circ_class_n(ipts,ivm,:)) |
|---|
| 422 | ENDIF |
|---|
| 423 | !- |
|---|
| 424 | |
|---|
| 425 | ELSEIF( (veget_max(ipts,ivm).LT. min_stomate) .AND. & |
|---|
| 426 | (SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .LT. min_stomate) .AND. & |
|---|
| 427 | (SUM(circ_class_n(ipts,ivm,:)) .LT. min_stomate) )THEN |
|---|
| 428 | |
|---|
| 429 | ! The PFT is not defined. Initialize the variables |
|---|
| 430 | ! to ensure that all goes well when moving and merging |
|---|
| 431 | ! different age classes. |
|---|
| 432 | circ_class_n(ipts,ivm,:) = zero |
|---|
| 433 | circ_class_biomass(ipts,ivm,:,:,icarbon) = zero |
|---|
| 434 | |
|---|
| 435 | ! Conditions are not fullfilled for establishing |
|---|
| 436 | ! new vegetation |
|---|
| 437 | establish = .FALSE. |
|---|
| 438 | recruite = .FALSE. |
|---|
| 439 | |
|---|
| 440 | ELSEIF (veget_max(ipts,ivm) .GT. min_stomate) THEN |
|---|
| 441 | |
|---|
| 442 | ! Conditions are not fullfilled for either establishing |
|---|
| 443 | ! a new young vegetation or having recruitment under an |
|---|
| 444 | ! exisisting vegetation of the same PFT. |
|---|
| 445 | establish = .FALSE. |
|---|
| 446 | recruite = .FALSE. |
|---|
| 447 | |
|---|
| 448 | ! Debug |
|---|
| 449 | IF (printlev_loc>=4) THEN |
|---|
| 450 | WRITE(numout,*) 'We are calling stomate_prescribe.f90 & |
|---|
| 451 | & but nothing should be done for PFT ',ivm |
|---|
| 452 | ENDIF |
|---|
| 453 | |
|---|
| 454 | ELSE |
|---|
| 455 | |
|---|
| 456 | ! Unexpected condition |
|---|
| 457 | CALL ipslerr_p (3,'stomate_prescribe.f90',& |
|---|
| 458 | 'Unexpected condition - time to rethink the', & |
|---|
| 459 | 'IF that choses between establish and recruite','') |
|---|
| 460 | |
|---|
| 461 | ENDIF |
|---|
| 462 | |
|---|
| 463 | !! 2.2 Calculate variables that will determine the size of the saplings |
|---|
| 464 | IF( establish .OR. recruite )THEN |
|---|
| 465 | |
|---|
| 466 | !! 2.2.1 Calculate c0_alloc |
|---|
| 467 | ! The calculation starts with c0_alloc, it only depends on PFT and |
|---|
| 468 | ! the effective longiveties. |
|---|
| 469 | c0_alloc(ivm) = calculate_c0_alloc(ipts, ivm, & |
|---|
| 470 | longevity_eff_root(ipts,ivm),longevity_eff_sap(ipts,ivm)) |
|---|
| 471 | |
|---|
| 472 | !! 2.2.2 Calculate k_latosa for newly established PFTs |
|---|
| 473 | ! If a new vegetation has to be established we recalculate |
|---|
| 474 | ! k_latosa from the parameter file (see note on physiological |
|---|
| 475 | ! memory of k_latosa below). If we add recruits we use the KF of |
|---|
| 476 | ! the existing vegetation. |
|---|
| 477 | IF (establish) THEN |
|---|
| 478 | |
|---|
| 479 | ! Initialize tree level variables |
|---|
| 480 | circ_class_n(ipts,ivm,:) = zero |
|---|
| 481 | circ_class_biomass(ipts,ivm,:,:,:) = zero |
|---|
| 482 | |
|---|
| 483 | ! Assume that there is no physiological memory for k_latosa |
|---|
| 484 | ! between different generations (this is probably true when |
|---|
| 485 | ! trees are planted but seems less likely for natural |
|---|
| 486 | ! regeneration. The first test that implemented memory has |
|---|
| 487 | ! caused problems with the LAI from the second generation |
|---|
| 488 | ! onwards. Nevertheless, we still use k_latosa_adapt. In case |
|---|
| 489 | ! someone want to try again in the future, the structure of |
|---|
| 490 | ! the code is ready to account for memory |
|---|
| 491 | k_latosa_adapt(ipts,ivm) = k_latosa_min(ivm) |
|---|
| 492 | |
|---|
| 493 | ! Debugging |
|---|
| 494 | IF (printlev_loc>=4) THEN |
|---|
| 495 | WRITE(numout,*) 'Establish a PFT (stomate_prescribe.f90):' |
|---|
| 496 | WRITE(numout,*) ' Imposing initial biomass for prescribed & |
|---|
| 497 | & trees, initial reserve mass for prescribed grasses.' |
|---|
| 498 | WRITE(numout,*) ' Declaring prescribed PFTs present.' |
|---|
| 499 | WRITE(numout,*) ' Grid point: ',ipts,' PFT Type: ',ivm |
|---|
| 500 | WRITE(numout,*) ' Establish ',establish,' Recruite ',recruite |
|---|
| 501 | ENDIF |
|---|
| 502 | !- |
|---|
| 503 | |
|---|
| 504 | ENDIF ! establish |
|---|
| 505 | |
|---|
| 506 | !! 2.2.3 Calculate Lightstress |
|---|
| 507 | ! Note that light and water stress have an opposite effect on KF. |
|---|
| 508 | ! Waterstress will decrease KF because more C should be allocated |
|---|
| 509 | ! to the roots. Light stress will increase KF because more C should |
|---|
| 510 | ! be allocated to the leaves. |
|---|
| 511 | ! Lightstress varies from 0 to 1 and is calculated from the canopy |
|---|
| 512 | ! structure (veget). |
|---|
| 513 | IF (establish) THEN |
|---|
| 514 | |
|---|
| 515 | ! Given that there is no vegetation at this point, the |
|---|
| 516 | ! lightstress cannot be calculated and is set to 1. |
|---|
| 517 | ! However, as soon as we grow a canopy it will experience |
|---|
| 518 | ! light stress and so KF should be adjusted. |
|---|
| 519 | lstress_fac(ipts,ivm) = un |
|---|
| 520 | |
|---|
| 521 | ELSEIF (recruite) THEN |
|---|
| 522 | |
|---|
| 523 | ! We already have vegetation, so we should do recruitment if |
|---|
| 524 | ! we want it for this PFT. First lets compute light stress based |
|---|
| 525 | ! on the fraction of light reaching the ground (level 1), |
|---|
| 526 | ! averaged over all daytime time steps (sinang>0) of the vegetative |
|---|
| 527 | ! period (gpp>0). This is based on Pgap and accounts for canopy |
|---|
| 528 | ! structure. |
|---|
| 529 | lstress_fac(ipts,ivm)=light_tran_to_floor_season(ipts,ivm) |
|---|
| 530 | |
|---|
| 531 | ! This value is accumulated in stomate.f90. It is accumulated |
|---|
| 532 | ! over an entire year and used only the last day of the year. |
|---|
| 533 | ! After it has been used, it can be reset. Note that |
|---|
| 534 | ! age_class_dist takes place before recruitment so it has be |
|---|
| 535 | ! moved from one PFT to another first. |
|---|
| 536 | light_tran_to_floor_season(ipts,ivm) = zero |
|---|
| 537 | |
|---|
| 538 | ! Debug |
|---|
| 539 | IF(printlev_loc>=4)THEN |
|---|
| 540 | WRITE(numout,*) 'lstress_fac in stomate_prescribe is ', & |
|---|
| 541 | lstress_fac(ipts,ivm)*100, ' (%)' |
|---|
| 542 | ENDIF |
|---|
| 543 | !- |
|---|
| 544 | |
|---|
| 545 | ENDIF ! establish/recruite for light |
|---|
| 546 | |
|---|
| 547 | ENDIF ! establish .OR. recruite for initialization |
|---|
| 548 | |
|---|
| 549 | !! 2.3 Initial vegetation has to be prescribed only when the vegetation is |
|---|
| 550 | ! static |
|---|
| 551 | IF (establish .OR. recruite) THEN |
|---|
| 552 | |
|---|
| 553 | !! 2.3.1 Initialize woody PFT's |
|---|
| 554 | ! Use veget_max to check whether the PFT is present. If the PFT |
|---|
| 555 | ! is present but it has no biomass, prescribe its biomass (gC m-{2}). |
|---|
| 556 | IF ( is_tree(ivm) ) THEN |
|---|
| 557 | |
|---|
| 558 | ! prescribe/recruitment branching for actual process |
|---|
| 559 | IF (establish) THEN |
|---|
| 560 | |
|---|
| 561 | ! Allocation factors |
|---|
| 562 | ! Sapwood to root ratio |
|---|
| 563 | ! Following Magnani et al. 2000 "In order to decreases hydraulic |
|---|
| 564 | ! resistance, the investment of carbon in fine roots or sapwood |
|---|
| 565 | ! yields to the plant very different returns, both because of |
|---|
| 566 | ! different hydraulic conductivities and because of the strong |
|---|
| 567 | ! impact of plant height on shoot resistance. On the other hand, |
|---|
| 568 | ! fine roots and sapwood have markedly different longevities and |
|---|
| 569 | ! the cost of production, discounted for turnover, will differ |
|---|
| 570 | ! accordingly. Optimal growth under hydraulic constraints requires |
|---|
| 571 | ! that the ratio of marginal hydraulic returns to marginal annual |
|---|
| 572 | ! cost for carbon investment in either roots or sapwood be the |
|---|
| 573 | ! same (Bloom, Chapin & Mooney 1985; Case & Fair 1989). This |
|---|
| 574 | ! is formalized in equation (13) and further derived to obtain |
|---|
| 575 | ! equation (17) in Magnani et al 2000. The latter is implemented |
|---|
| 576 | ! here.Pipe_density is given in gC/m-3, convert to kg/m-3. And |
|---|
| 577 | ! apply equation (17) in Magnani et al 2000. Note that c0_alloc |
|---|
| 578 | ! was calculated at the start of this routine. The calculation |
|---|
| 579 | ! itself is done in function_library |
|---|
| 580 | |
|---|
| 581 | ! Calculate leaf area to sapwood area |
|---|
| 582 | ! To be consistent with the hydraulic limitations and pipe theory, |
|---|
| 583 | ! k_latosa is calculated from equation (18) in Magnani et al. |
|---|
| 584 | ! To do so, total hydraulic resistance and tree height need to |
|---|
| 585 | ! known. This poses a problem as the resistance depends on the leaf |
|---|
| 586 | ! area and the leaf area on the resistance. There is no independent |
|---|
| 587 | ! equation and equations 12 and 18 depend on each other and |
|---|
| 588 | ! substitution would be circular. Hence prescribed k_latosa values |
|---|
| 589 | ! were obtained from observational records and are given in |
|---|
| 590 | ! mtc_parameters.f90. |
|---|
| 591 | |
|---|
| 592 | ! The relationship between height and k_latosa as reported in |
|---|
| 593 | ! McDowell et al 2002 and Novick et al 2009 could be implemented |
|---|
| 594 | ! to adjust k_latosa for the height of the stand. This did NOT |
|---|
| 595 | ! result in a realistic model behavior |
|---|
| 596 | ! k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * & |
|---|
| 597 | ! (k_latosa_max(ivm) - (latosa_height(ivm) * & |
|---|
| 598 | ! (SUM( nb_trees_i(:) * ave_tree_height(:) ) / & |
|---|
| 599 | ! SUM( nb_trees_i(:) )))) |
|---|
| 600 | |
|---|
| 601 | ! Also k_latosa has been reported to be a function of CO2 |
|---|
| 602 | ! concentration (Atwell et al. 2003, Tree Physiology, 23:13-21 |
|---|
| 603 | ! and Pakati et al. 2000, Global Change Biology, 6:889-897). |
|---|
| 604 | ! This effect is not accounted for in the current code |
|---|
| 605 | |
|---|
| 606 | ! Finally, k_latosa is also reported to be a function of |
|---|
| 607 | ! diameter (i.e. stand thinning, Simonin et al 2006, Tree |
|---|
| 608 | ! Physiology, 26:493-503). Here the relationship with thinning was |
|---|
| 609 | ! interpreted as a realtionship with light stress. Note that light |
|---|
| 610 | ! stress cannot be calculated at this time in the model because |
|---|
| 611 | ! there is no canopy (that's why we are in prescribe!) and |
|---|
| 612 | ! so there is no lightstress. lstress was therefore set to one |
|---|
| 613 | ! (see above). We prefered to keep this redundancy in the code |
|---|
| 614 | ! because it makes it clear that k_latosa is calculated in the |
|---|
| 615 | ! same way in prescribe.f90 and growth_fun_all.f90 |
|---|
| 616 | ! k_latosa(ipts,ivm) = (k_latosa_adapt(ipts,ivm) + & |
|---|
| 617 | ! (lstress_fac(ipts,ivm) * & |
|---|
| 618 | ! (k_latosa_max(ivm)-k_latosa_min(ivm)))) |
|---|
| 619 | |
|---|
| 620 | ! +++CHECK+++ |
|---|
| 621 | ! Ideally waterstress and lightstress should both |
|---|
| 622 | ! affect k_latosa. One of the following approaches could |
|---|
| 623 | ! be used as a starting point. |
|---|
| 624 | ! k_latosa(ipts,ivm) = k_latosa_min(ivm) + & |
|---|
| 625 | ! (wstress_fac(ipts,ivm) * lstress_fac(ipts,ivm) * & |
|---|
| 626 | ! (k_latosa_max(ivm)-k_latosa_min(ivm))) |
|---|
| 627 | ! k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * (k_latosa_min(ivm) + & |
|---|
| 628 | ! (lstress_fac(ipts,ivm) * & |
|---|
| 629 | ! (k_latosa_max(ivm)-k_latosa_min(ivm)))) |
|---|
| 630 | ! ++++++++++++ |
|---|
| 631 | |
|---|
| 632 | ! Calculate conversion coefficient for sapwood area to leaf |
|---|
| 633 | ! area |
|---|
| 634 | ! (1) The scaling parameter between leaf and sapwood mass |
|---|
| 635 | ! is derived from LA_ind = k_latosa * SA_ind, where |
|---|
| 636 | ! LA_ind = leaf area of an individual, SA_ind is the sapwood |
|---|
| 637 | ! area of an individual and k_latosa a pipe-model parameter |
|---|
| 638 | ! (2) LA_ind = Cl * sla |
|---|
| 639 | ! (3) Cs = SA_ind * height * wooddensity * tree_ff |
|---|
| 640 | ! Substitute (2) and (3) in (1) |
|---|
| 641 | ! Cl = Cs * k1 / (wooddensity * sla * tree_ff * height) |
|---|
| 642 | ! Cl = Cs*KF/height, where KF is in (m) |
|---|
| 643 | ! KF is passed to the allocation routine and it is saved in the |
|---|
| 644 | ! restart file. |
|---|
| 645 | !+++CHECK+++ |
|---|
| 646 | ! At one point it looked like a good idea to take the max of two |
|---|
| 647 | ! options but by doing so we cannot recalculate KF in phenology. |
|---|
| 648 | ! It has not been confirmed that this is really a problem. Use the |
|---|
| 649 | ! simpelest approach but leave the alternative in the code as a |
|---|
| 650 | ! suggestion of a possible solution in case something goes wrong. |
|---|
| 651 | !!$ k_latosa_tmp = MAX(k_latosa_min(ivm),k_latosa_adapt(ipts,ivm) + & |
|---|
| 652 | !!$ (lstress_fac(ipts,ivm) * & |
|---|
| 653 | !!$ (k_latosa_max(ivm)-k_latosa_min(ivm)))) |
|---|
| 654 | !+++++++++++ |
|---|
| 655 | k_latosa_tmp = k_latosa_adapt(ipts,ivm) + & |
|---|
| 656 | (lstress_fac(ipts,ivm) * & |
|---|
| 657 | (k_latosa_max(ivm)-k_latosa_min(ivm))) |
|---|
| 658 | |
|---|
| 659 | IF (sla_dyn) THEN |
|---|
| 660 | KF(ipts,ivm) = k_latosa_tmp / & |
|---|
| 661 | (slainit(ivm) * pipe_density(ivm) * tree_ff(ivm)) |
|---|
| 662 | ELSE |
|---|
| 663 | KF(ipts,ivm) = k_latosa_tmp / & |
|---|
| 664 | (sla(ivm) * pipe_density(ivm) * tree_ff(ivm)) |
|---|
| 665 | END IF |
|---|
| 666 | |
|---|
| 667 | ! Multiple approaches for initializing the forest have been tried. |
|---|
| 668 | ! 1) Using self-thinning solely. This method plants an unrealistically high |
|---|
| 669 | ! amount of individuals which results in too little growth per individual |
|---|
| 670 | ! at the start. Self-thinning lines are derived from forest inventory data |
|---|
| 671 | ! which has size threshold of trees for data collection, so we should not |
|---|
| 672 | ! rely on the tails of the self-thinning relationships which depend on |
|---|
| 673 | ! observations of saplings. Other approaches were therefore necessary to |
|---|
| 674 | ! better initialize the model. |
|---|
| 675 | ! 2) Constant mortality (as in r2566) until the trees are big enough to |
|---|
| 676 | ! reach the self-thinning line. The initial number of individuals (IND) was |
|---|
| 677 | ! prescribed per PFT and the background mortality was prescribed across all |
|---|
| 678 | ! PFTs. It worked mostly OK in the previous version, but some problems (very |
|---|
| 679 | ! slow growth at the start) was observed for isolated pixels/species. There is |
|---|
| 680 | ! no control over when self-thinning is applied. This approach overlooks the |
|---|
| 681 | ! fact that there are certain sections of the self-thinning relationship that |
|---|
| 682 | ! we trust better than other sections. |
|---|
| 683 | ! 3) Control by initial diameter and RDI. The initial IND was calculated |
|---|
| 684 | ! using prescribed RDI. It worked well with most of species, but some species |
|---|
| 685 | ! showed problems because the initial number of individuals was already in |
|---|
| 686 | ! the tail of the self-thinning relationship which resulted in very slow |
|---|
| 687 | ! initial growth. For other species, the initial diameter was quite high |
|---|
| 688 | ! (> 5cm) for the maximum number of trees (set to 14,000) which implied that |
|---|
| 689 | ! we were planting 10+ years old stands. It also made the first age class |
|---|
| 690 | ! almost redundant. This increased the risk of unrealistic carbon balances |
|---|
| 691 | ! (too fast growing because year 1 in the model contained 10 yearsâ worth of |
|---|
| 692 | ! carbon) and unrealistic biophysics because the initial LAI and tree height |
|---|
| 693 | ! were not representing a stand of saplings. |
|---|
| 694 | ! 4) Calculate initial IND based on fixed initial diameter and LAI and IND |
|---|
| 695 | ! decreases linearly towards the first trustworthy point of the self-thinning |
|---|
| 696 | ! relationship. It is a similar idea with constant mortality but we could control |
|---|
| 697 | ! the number of individuals and their diameter for which the model switches to |
|---|
| 698 | ! self-thinning. The whole IND~D line is composed of self-thinning line and |
|---|
| 699 | ! the linear line connecting the prescribed threshold of IND we think trustable |
|---|
| 700 | ! based on a dataset and initial N calculated from diameter and LAI. The |
|---|
| 701 | ! problem was it has different gradients by species due to self-thinning line |
|---|
| 702 | ! and parameters. Species with low gradient didn't grow well because too many |
|---|
| 703 | ! individuals slowed down the growth per individual. |
|---|
| 704 | ! 5) Calculate initial IND based on fixed gradient and a trustworthy starting |
|---|
| 705 | ! point on the self-thinning relationship. A thinning line has the same shape |
|---|
| 706 | ! with 4) but the gradient of linear line is fixed in this case. It worked |
|---|
| 707 | ! well with all species. |
|---|
| 708 | ! 6) For managed forest a yield based RDI approch is used to |
|---|
| 709 | ! thin the trees at low density in order to fasten the young forest stage. |
|---|
| 710 | ! !!WARNING!! if the RDI upper target is too low the forest cannot accumalate |
|---|
| 711 | ! enough sugar for the next budburst due to a low LAI. To avoid that we decide |
|---|
| 712 | ! to set a thresold below which a fix value of 0.4 for RDI upper target is used. |
|---|
| 713 | IF(forest_managed(ipts,ivm) .EQ. ifm_none) THEN |
|---|
| 714 | !Approach 6 |
|---|
| 715 | n_init(ivm) = & |
|---|
| 716 | MIN(Nmax(qmd_init(ivm)*m_to_cm,ivm)* & |
|---|
| 717 | MAX(MIN(a_rdi_upper_unman(ivm)+ qmd_init(ivm)* & |
|---|
| 718 | b_rdi_upper_unman(ivm),c_rdi_upper_unman(ivm)),& |
|---|
| 719 | d_rdi_upper_unman(ivm)),nmaxplants(ivm)) |
|---|
| 720 | |
|---|
| 721 | ELSE |
|---|
| 722 | !Approach 6 |
|---|
| 723 | n_init(ivm) = & |
|---|
| 724 | MIN(Nmax(qmd_init(ivm)*m_to_cm,ivm)* & |
|---|
| 725 | MAX(MIN(a_rdi_upper_man(ivm)+ qmd_init(ivm)* & |
|---|
| 726 | b_rdi_upper_man(ivm),c_rdi_upper_man(ivm)),& |
|---|
| 727 | d_rdi_upper_man(ivm)),nmaxplants(ivm)) |
|---|
| 728 | ENDIF |
|---|
| 729 | n_init(ivm) = n_init(ivm)*ha_to_m2 |
|---|
| 730 | |
|---|
| 731 | ! Debug |
|---|
| 732 | IF (printlev_loc.GT.4) THEN |
|---|
| 733 | WRITE(numout,*) 'n_init (trees/m2), ', ivm, n_init(ivm) |
|---|
| 734 | END IF |
|---|
| 735 | !- |
|---|
| 736 | |
|---|
| 737 | DO icir = 1,ncirc |
|---|
| 738 | ! Make use of dia_init (calculated in stomate.f90 as a parameter) |
|---|
| 739 | ! do calculate the height and wood mass of the saplings of a newly |
|---|
| 740 | ! planted PFT. |
|---|
| 741 | ave_tree_height(icir) = pipe_tune2(ivm)* & |
|---|
| 742 | (dia_init(ivm,icir)**pipe_tune3(ivm)) |
|---|
| 743 | wood_init = ( ave_tree_height(icir) * pi / 4. * & |
|---|
| 744 | (dia_init(ivm,icir)) ** 2. ) * pipe_density(ivm) * tree_ff(ivm) |
|---|
| 745 | |
|---|
| 746 | CALL make_saplings(bm_sapl_2D(ipts,ivm,icir,:,:), & |
|---|
| 747 | dia_init(ivm,icir), icir, ivm, & |
|---|
| 748 | ipts, KF(ipts,ivm), c0_alloc(ivm),plant_status(ipts,ivm), & |
|---|
| 749 | establish, cn_leaf_init_2D(ipts,ivm),bypass_lm(icir)) |
|---|
| 750 | ENDDO |
|---|
| 751 | |
|---|
| 752 | ! Prescribe the newly planted population. |
|---|
| 753 | circ_class_biomass(ipts,ivm,:,:,:) = bm_sapl_2D(ipts,ivm,:,:,:) |
|---|
| 754 | circ_class_n(ipts,ivm,:) = circ_class_dist(:) / & |
|---|
| 755 | SUM(circ_class_dist(:)) * n_init(ivm) |
|---|
| 756 | |
|---|
| 757 | ! Debug |
|---|
| 758 | IF (printlev_loc>=4) THEN |
|---|
| 759 | DO icir=1,ncirc |
|---|
| 760 | WRITE(numout,*)'icir, circ_class_n, ', icir, & |
|---|
| 761 | circ_class_n(ipts,ivm,icir) |
|---|
| 762 | WRITE(numout,*)'Initial biomass distribution - carbon' |
|---|
| 763 | WRITE(numout,*) circ_class_biomass(ipts,ivm,icir,:,icarbon) |
|---|
| 764 | WRITE(numout,*) circ_class_biomass(ipts,ivm,icir,:,initrogen) |
|---|
| 765 | ENDDO |
|---|
| 766 | WRITE(numout,*)' qmd_init', qmd_init(ivm) |
|---|
| 767 | WRITE(numout,*)' n_init(ivm)', n_init(ivm) |
|---|
| 768 | WRITE(numout,*)' Maximun number of saplings', Nmax(qmd_init(ivm),ivm) |
|---|
| 769 | END IF |
|---|
| 770 | !- |
|---|
| 771 | |
|---|
| 772 | ELSEIF (recruite) THEN |
|---|
| 773 | |
|---|
| 774 | !! Calculate the number of recruits |
|---|
| 775 | ! Recruitment is based on the light stress factor [0 1] computed |
|---|
| 776 | ! above. It makes use of the functional response of recruitment |
|---|
| 777 | ! to light described by Ruger et al (2009, J. of Ecol., |
|---|
| 778 | ! doi: 10.1111/j.1365-2745.2009.01552.x) in the 50-ha plot of Barro |
|---|
| 779 | ! Colorado Island (BCI) in Panama: |
|---|
| 780 | ! Log10(Nrecruits)= A + B * ( Log10(lstress_fac) - MU ) |
|---|
| 781 | ! |
|---|
| 782 | ! Where Nrecruits is the number of recruits for an area of 25 m2. |
|---|
| 783 | ! The parameter A is the intercept (mean log10 recruits expected), |
|---|
| 784 | ! B the slope (functional response), and MU=Log10(0.02)=-1.7 |
|---|
| 785 | ! Here the recruitment response to light can be negative (B<0), |
|---|
| 786 | ! decelerating (0<B<1), accelerating (B>1) or linear (B=1). |
|---|
| 787 | ! MU is needed here to express the fraction of light between 0 and 1 |
|---|
| 788 | ! centred on mean light 2% at BCI. The linearized (log-log) power |
|---|
| 789 | ! model is implemented here for each PFT following the original setup |
|---|
| 790 | ! of Ruger et al 2009 to ensure realistic values. IT COULD EASILY BE |
|---|
| 791 | ! SIMPLIFIED. The predicted number of recruits is for an area of |
|---|
| 792 | ! 25 m2 so we have to compute the antilogarithm and divide by 25 to |
|---|
| 793 | ! get the estimate by 1 m2 which is the unit used in ORCHIDEE. |
|---|
| 794 | ! |
|---|
| 795 | ! Add min_stomate (a very small number) to lstress_fac to avoid numerical |
|---|
| 796 | ! issues with LOG10(zero) being undefined. |
|---|
| 797 | IF (recruitment_alpha(ivm).LT.min_stomate .AND. & |
|---|
| 798 | recruitment_beta(ivm).LT.min_stomate) THEN |
|---|
| 799 | |
|---|
| 800 | ! Sanity check. If the users sets both values to zero |
|---|
| 801 | ! new_ind is set to zero. The simulation should be |
|---|
| 802 | ! identical to a simulation without recruitment. |
|---|
| 803 | new_ind(ipts,ivm) = zero |
|---|
| 804 | |
|---|
| 805 | ELSE |
|---|
| 806 | |
|---|
| 807 | ! Calculate the number of recruits m2 day-1 |
|---|
| 808 | new_ind(ipts,ivm) = ( 10.**( recruitment_alpha(ivm) + & |
|---|
| 809 | recruitment_beta(ivm) * & |
|---|
| 810 | (LOG10(lstress_fac(ipts,ivm)+min_stomate) - & |
|---|
| 811 | LOG10(0.02)) ) )/25. |
|---|
| 812 | |
|---|
| 813 | ENDIF |
|---|
| 814 | |
|---|
| 815 | ! Update the number of individuals |
|---|
| 816 | tmp_ind = SUM(circ_class_n(ipts,ivm,:)) |
|---|
| 817 | old_ind = tmp_ind |
|---|
| 818 | tmp_ind = tmp_ind + new_ind(ipts,ivm) |
|---|
| 819 | |
|---|
| 820 | ! Debug |
|---|
| 821 | IF(printlev_loc>=4)THEN |
|---|
| 822 | WRITE(numout,*) 'PFT, ',ivm |
|---|
| 823 | WRITE(numout,*) 'canopy cover (LAI), ',& |
|---|
| 824 | cc_to_lai(circ_class_biomass(ipts,ivm,:,ileaf,icarbon),& |
|---|
| 825 | circ_class_n(ipts,ivm,:),ivm) |
|---|
| 826 | WRITE(numout,*) 'used canopy_gap (%) is, ', lstress_fac(ipts,ivm)*100 |
|---|
| 827 | WRITE(numout,*) 'Original ind is, ', (tmp_ind - new_ind(ipts,ivm))*m2_to_ha |
|---|
| 828 | WRITE(numout,*) 'new_ind is, ', new_ind(ipts,ivm)*m2_to_ha |
|---|
| 829 | WRITE(numout,*) 'ind+new_ind is, ', tmp_ind*m2_to_ha |
|---|
| 830 | WRITE(numout,*) 'Dg (m) ', qmd_init(ivm) |
|---|
| 831 | ind_init = Nmax(qmd_init(ivm)*m_to_cm,ivm) |
|---|
| 832 | WRITE(numout,*) 'Nmax, ',ind_init |
|---|
| 833 | ENDIF |
|---|
| 834 | !- |
|---|
| 835 | |
|---|
| 836 | !! Calculate the dimensions of the recruits |
|---|
| 837 | ! The height of the saplings is prescribed and determines the reserves |
|---|
| 838 | ! which are especially important for deciduous species, which need to |
|---|
| 839 | ! survive on their reserves for the first couple of months to a whole |
|---|
| 840 | ! year because the phenology scheme requires annual mean values to get |
|---|
| 841 | ! started. Calculate sapling diameter and biomass for prescribed sapling |
|---|
| 842 | ! height (above). Note that the dimensions of a newly planted sapling and |
|---|
| 843 | ! a recruit are not the same. |
|---|
| 844 | dia_init_recr = ( recruitment_height(ivm) / pipe_tune2(ivm) ) ** & |
|---|
| 845 | ( 1. / pipe_tune3(ivm) ) |
|---|
| 846 | wood_init = (recruitment_height(ivm) * pi / 4. * (dia_init_recr) ** 2. ) * & |
|---|
| 847 | pipe_density(ivm) * tree_ff(ivm) |
|---|
| 848 | |
|---|
| 849 | ! Use the variables from the pipe model and the allometric |
|---|
| 850 | ! relationships to calculate the different biomass components |
|---|
| 851 | ! of the saplings in each diameter class. |
|---|
| 852 | CALL make_saplings(bm_sapl_2D(ipts,ivm,1,:,:), dia_init_recr, & |
|---|
| 853 | 1, ivm, ipts, KF(ipts,ivm), c0_alloc(ivm), plant_status(ipts,ivm), & |
|---|
| 854 | establish, cn_leaf_init_2D(ipts,ivm),bypass_lm(1)) |
|---|
| 855 | |
|---|
| 856 | ! Debug |
|---|
| 857 | IF (printlev_loc>=4) THEN |
|---|
| 858 | WRITE(numout,*)'Prescribed recruitment height (m): ',recruitment_height(ivm) |
|---|
| 859 | WRITE(numout,*)'Recruitment diameter (m) for prescribed height: ', dia_init |
|---|
| 860 | WRITE(numout,*)'Wood_init recruitment: ',wood_init |
|---|
| 861 | WRITE(numout,*)'ind_init ',ind_init |
|---|
| 862 | WRITE(numout,*)'bm_sapl after calling make_saplings: ',bm_sapl_2D(ipts,ivm,1,:,:) |
|---|
| 863 | ENDIF |
|---|
| 864 | !- |
|---|
| 865 | |
|---|
| 866 | ! Update the number of individuals in the first size class. All |
|---|
| 867 | ! recruitment occurs in the circ_class with the smallest diameter |
|---|
| 868 | old_circ_class_n(ipts,ivm,:) = circ_class_n(ipts,ivm,:) |
|---|
| 869 | circ_class_n(ipts,ivm,1) = circ_class_n(ipts,ivm,1) + new_ind(ipts,ivm) |
|---|
| 870 | |
|---|
| 871 | !! Determine the biomass in the first circumference class. |
|---|
| 872 | ! This is a dillution approach were an average tree is made by making use |
|---|
| 873 | ! of the biomass in the available trees and the biomass in the recruits. |
|---|
| 874 | ! This is a simple mean weighted by the number of trees. The model needs |
|---|
| 875 | ! the biomass in an individual tree (gC tree-1) |
|---|
| 876 | DO iele = 1,nelements |
|---|
| 877 | old_circ_class_biomass(ipts,ivm,:,:,iele) = & |
|---|
| 878 | circ_class_biomass(ipts,ivm,:,:,iele) |
|---|
| 879 | circ_class_biomass(ipts,ivm,1,:,iele) = & |
|---|
| 880 | (circ_class_biomass(ipts,ivm,1,:,iele) * & |
|---|
| 881 | old_circ_class_n(ipts,ivm,1) + & |
|---|
| 882 | bm_sapl_2D(ipts,ivm,1,:,iele) * new_ind(ipts,ivm) ) / & |
|---|
| 883 | circ_class_n(ipts,ivm,1) |
|---|
| 884 | ENDDO |
|---|
| 885 | |
|---|
| 886 | ! Debug |
|---|
| 887 | IF (printlev_loc>=4) THEN |
|---|
| 888 | WRITE(numout,*) "Canopy gap (%)", lstress_fac(ipts,ivm)*100 |
|---|
| 889 | WRITE(numout,*) "Original density (N/ha): ",& |
|---|
| 890 | old_circ_class_n(ipts,ivm,:)*m2_to_ha |
|---|
| 891 | WRITE(numout,*) "New density (N/ha): ",circ_class_n(ipts,ivm,:)*m2_to_ha |
|---|
| 892 | WRITE(numout,*) "Recruited saplings per ha in first size class ", & |
|---|
| 893 | new_ind(ipts,ivm) * m2_to_ha |
|---|
| 894 | WRITE(numout,*) "Original biomass per individual", & |
|---|
| 895 | SUM(old_circ_class_biomass(ipts,ivm,:,:,icarbon),2),& |
|---|
| 896 | SUM(old_circ_class_biomass(ipts,ivm,:,:,initrogen),2) |
|---|
| 897 | WRITE(numout,*) "New biomass per individual", & |
|---|
| 898 | SUM( circ_class_biomass(ipts,ivm,:,:,icarbon),2 ),& |
|---|
| 899 | SUM( circ_class_biomass(ipts,ivm,:,:,initrogen),2 ) |
|---|
| 900 | WRITE(numout,*) "Biomass of new saplings gC/sapling", & |
|---|
| 901 | bm_sapl_2D(ipts,ivm,1,:,icarbon) |
|---|
| 902 | WRITE(numout,*) "Biomass of new saplings gN/sapling", & |
|---|
| 903 | bm_sapl_2D(ipts,ivm,1,:,initrogen) |
|---|
| 904 | WRITE(numout,*) "KF ", KF(ipts,ivm) |
|---|
| 905 | ENDIF |
|---|
| 906 | !- |
|---|
| 907 | |
|---|
| 908 | ENDIF ! end prescribe/recruitment block |
|---|
| 909 | |
|---|
| 910 | IF (establish) THEN |
|---|
| 911 | |
|---|
| 912 | ! Set leaf age classes, all leaves are current year leaves |
|---|
| 913 | leaf_frac(ipts,ivm,:) = zero |
|---|
| 914 | leaf_frac(ipts,ivm,1) = un |
|---|
| 915 | |
|---|
| 916 | ! Set time since last beginning of growing season but only |
|---|
| 917 | ! for the first day of the whole simulation. When the model |
|---|
| 918 | ! is initialized when_growthinit is set to undef. In subsequent |
|---|
| 919 | ! time steps it should have a value. For trees without phenology |
|---|
| 920 | ! the growing season starts at the moment the PFT is prescribed |
|---|
| 921 | IF (when_growthinit(ipts,ivm) .EQ. undef) THEN |
|---|
| 922 | when_growthinit(ipts,ivm) = 200 |
|---|
| 923 | ENDIF |
|---|
| 924 | |
|---|
| 925 | ENDIF |
|---|
| 926 | |
|---|
| 927 | ! The carbon and nitrogen to build the saplings is taken from the |
|---|
| 928 | ! atmosphere, keep track of amount to calculate the C-balance |
|---|
| 929 | ! closure |
|---|
| 930 | IF (establish) THEN |
|---|
| 931 | |
|---|
| 932 | ! The whole biomass of this PFT was taken from the |
|---|
| 933 | ! atmosphere |
|---|
| 934 | atm_to_bm(ipts,ivm,:) = atm_to_bm(ipts,ivm,:) + & |
|---|
| 935 | (SUM(cc_to_biomass(npts, ivm, & |
|---|
| 936 | circ_class_biomass(ipts,ivm,:,:,:), & |
|---|
| 937 | circ_class_n(ipts,ivm,:)),1) / dt ) |
|---|
| 938 | |
|---|
| 939 | ELSEIF (recruite) THEN |
|---|
| 940 | |
|---|
| 941 | ! Only the increase in biomass of this PFT was taken |
|---|
| 942 | ! from the atmosphere |
|---|
| 943 | atm_to_bm(ipts,ivm,:) = atm_to_bm(ipts,ivm,:) + & |
|---|
| 944 | (SUM(cc_to_biomass(npts, ivm, & |
|---|
| 945 | circ_class_biomass(ipts,ivm,:,:,:), & |
|---|
| 946 | circ_class_n(ipts,ivm,:)) - & |
|---|
| 947 | old_biomass(ipts,ivm,:,:),1) / dt ) |
|---|
| 948 | |
|---|
| 949 | ENDIF |
|---|
| 950 | |
|---|
| 951 | ! Check whether the order was preserved. We only want to |
|---|
| 952 | ! preserve the order for the carbon biomass. Nitrogen follows |
|---|
| 953 | ! carbon. We expect that the biomass of an individual tree |
|---|
| 954 | ! increases with an increasing circ class. This is probably |
|---|
| 955 | ! not important for the correct functioning of the model. So, |
|---|
| 956 | ! if this turns out to be computationally expensive it is |
|---|
| 957 | ! worth trying without. We try to preserve the order as |
|---|
| 958 | ! an additional mean to control the quality of the simulations. |
|---|
| 959 | ! A strict order also helps to interpret the output. Hence, |
|---|
| 960 | ! if the order was not preserved we will sort the biomasses. |
|---|
| 961 | ! Note that this subroutine first checks whether the order |
|---|
| 962 | ! was preseverd. |
|---|
| 963 | circ_class_biomass_new = circ_class_biomass(ipts,ivm,:,:,:) |
|---|
| 964 | circ_class_n_new = circ_class_n(ipts,ivm,:) |
|---|
| 965 | CALL sort_circ_class_biomass(circ_class_biomass_new, & |
|---|
| 966 | circ_class_n_new) |
|---|
| 967 | circ_class_biomass(ipts,ivm,:,:,:) = circ_class_biomass_new |
|---|
| 968 | circ_class_n(ipts,ivm,:) = circ_class_n_new |
|---|
| 969 | |
|---|
| 970 | ENDIF ! tree(ivm) |
|---|
| 971 | |
|---|
| 972 | !! 2.3.2 Initialize grassy PFTs |
|---|
| 973 | ! Use veget_max to check whether the PFT is present. If the PFT is |
|---|
| 974 | ! present but it has no biomass, prescribe its biomass (gC m-{2}). |
|---|
| 975 | ! It is assumed that at day 1 grasses have all their biomass in the |
|---|
| 976 | ! reserve pool. The criteria exclude crops. Crops are no longer |
|---|
| 977 | ! prescribed but planted the day that begin_leaves is true. |
|---|
| 978 | IF (establish) THEN |
|---|
| 979 | |
|---|
| 980 | IF ( ( .NOT. is_tree(ivm) ) .AND. & |
|---|
| 981 | ( natural(ivm) ) .AND. & |
|---|
| 982 | ( veget_max(ipts,ivm) .GT. min_stomate ) .AND. & |
|---|
| 983 | ( SUM(circ_class_biomass(ipts,ivm,1,:,icarbon) ) .LE. min_stomate ) ) THEN |
|---|
| 984 | |
|---|
| 985 | ! Debug |
|---|
| 986 | IF(printlev_loc>=4)THEN |
|---|
| 987 | WRITE(numout,*) 'We will prescribe a new grass '//& |
|---|
| 988 | 'vegetation, the old one died' |
|---|
| 989 | ENDIF |
|---|
| 990 | ! - |
|---|
| 991 | |
|---|
| 992 | ! For grasses we assume that an individual grass is 1 m2 of grass. |
|---|
| 993 | ! This is set in nmaxplants in pft_parameters.f90. It could be set |
|---|
| 994 | ! to another value but with the current code this should not have |
|---|
| 995 | ! any meaning. The grassland does not necessarily covers the whole |
|---|
| 996 | ! 1 m2 so adjust for the canopy cover |
|---|
| 997 | tmp_ind = nmaxplants(ivm) * ha_to_m2 * canopy_cover(ivm) |
|---|
| 998 | |
|---|
| 999 | !+++CHECK+++ |
|---|
| 1000 | ! We do not deal with circumference classes for grasses |
|---|
| 1001 | ! and crops, but we want to keep the arrays the same as for the trees |
|---|
| 1002 | ! so we put the young grass information into the first circumference |
|---|
| 1003 | ! class. Calculate the sapwood to leaf mass in a similar way as has |
|---|
| 1004 | ! been done for trees. For trees this approach had been justified by |
|---|
| 1005 | ! observations. For grasses such justification is not supported by |
|---|
| 1006 | ! observations but we didn't try to find it. Needs more work by |
|---|
| 1007 | ! someone interested in grasses. There might be a more elegant |
|---|
| 1008 | ! solution making use of a well observed parameter. |
|---|
| 1009 | ! The mass of the structural carbon relates to the mass of the leaves |
|---|
| 1010 | ! through a prescribed parameter ::k_latosa |
|---|
| 1011 | !+++CHECK+++ |
|---|
| 1012 | ! At one point it looked like a good idea to take the max of two |
|---|
| 1013 | ! options but by doing so we cannot recalculate KF in phenology. |
|---|
| 1014 | ! It has not been confirmed that this is really a problem. Use the |
|---|
| 1015 | ! simpelest approach but leave the alternative in the code as a |
|---|
| 1016 | ! suggestion of a possible solution in case something goes wrong. |
|---|
| 1017 | !!$ k_latosa_tmp = MAX(k_latosa_min(ivm),k_latosa_adapt(ipts,ivm) + & |
|---|
| 1018 | !!$ (lstress_fac(ipts,ivm) * & |
|---|
| 1019 | !!$ (k_latosa_max(ivm)-k_latosa_min(ivm)))) |
|---|
| 1020 | !+++++++++++ |
|---|
| 1021 | k_latosa_tmp = k_latosa_adapt(ipts,ivm) + lstress_fac(ipts,ivm) * & |
|---|
| 1022 | (k_latosa_max(ivm)-k_latosa_min(ivm)) |
|---|
| 1023 | |
|---|
| 1024 | IF (sla_dyn) THEN |
|---|
| 1025 | KF(ipts,ivm) = k_latosa_tmp / & |
|---|
| 1026 | (slainit(ivm) * tree_ff(ivm) * pipe_density(ivm)) |
|---|
| 1027 | ELSE |
|---|
| 1028 | KF(ipts,ivm) = k_latosa_tmp / & |
|---|
| 1029 | (sla(ivm) * tree_ff(ivm) * pipe_density(ivm)) |
|---|
| 1030 | END IF |
|---|
| 1031 | |
|---|
| 1032 | ! Calculate leaf to root area |
|---|
| 1033 | LF(ipts,ivm) = c0_alloc(ivm) * KF(ipts,ivm) |
|---|
| 1034 | |
|---|
| 1035 | ! Debug |
|---|
| 1036 | IF(printlev_loc>=4.AND. test_pft == ivm)THEN |
|---|
| 1037 | WRITE(numout,*) 'KF, ', KF(ipts,ivm) |
|---|
| 1038 | WRITE(numout,*) 'LF, c0_alloc, ', c0_alloc(ivm) * KF(ipts,ivm), & |
|---|
| 1039 | c0_alloc(ivm) |
|---|
| 1040 | ENDIF |
|---|
| 1041 | !- |
|---|
| 1042 | |
|---|
| 1043 | ! initialize everything to make sure there are not random values |
|---|
| 1044 | ! floating around for ncirc = 1 |
|---|
| 1045 | bm_sapl_2D(ipts,ivm,1,:,:) = val_exp |
|---|
| 1046 | |
|---|
| 1047 | ! Similar as for trees, the initial height of the vegetation was |
|---|
| 1048 | ! defined. By using lai_to_biomass dynamic sla is accounted for. |
|---|
| 1049 | bm_sapl_2D(ipts,ivm,1,ileaf,icarbon) = lai_to_biomass(& |
|---|
| 1050 | height_init(ivm)/lai_to_height(ivm),ivm) |
|---|
| 1051 | |
|---|
| 1052 | ! Use allometric relationships to define the root mass based on |
|---|
| 1053 | ! leaf mass. Some 'sapwood mass' is needed to store the reserves. |
|---|
| 1054 | ! An arbitrairy fraction of 5% was used |
|---|
| 1055 | bm_sapl_2D(ipts,ivm,1,iroot,icarbon) = bm_sapl_2D(ipts,ivm,1,ileaf,icarbon) / & |
|---|
| 1056 | LF(ipts,ivm) |
|---|
| 1057 | bm_sapl_2D(ipts,ivm,1,isapabove,icarbon) = bm_sapl_2D(ipts,ivm,1,ileaf,icarbon) / & |
|---|
| 1058 | KF(ipts,ivm) |
|---|
| 1059 | |
|---|
| 1060 | ! Some of the biomass components that exist for trees are undefined |
|---|
| 1061 | ! for grasses |
|---|
| 1062 | bm_sapl_2D(ipts,ivm,1,isapbelow,icarbon) = zero |
|---|
| 1063 | bm_sapl_2D(ipts,ivm,1,iheartabove,icarbon) = zero |
|---|
| 1064 | bm_sapl_2D(ipts,ivm,1,iheartbelow,icarbon) = zero |
|---|
| 1065 | bm_sapl_2D(ipts,ivm,1,ifruit,icarbon) = zero |
|---|
| 1066 | bm_sapl_2D(ipts,ivm,1,icarbres,icarbon) = zero |
|---|
| 1067 | |
|---|
| 1068 | ! Pools that are defined in the same way for trees and grasses |
|---|
| 1069 | bm_sapl_2D(ipts,ivm,1,ilabile,icarbon) = labile_to_total * & |
|---|
| 1070 | (bm_sapl_2D(ipts,ivm,1,ileaf,icarbon) + & |
|---|
| 1071 | fcn_root(ivm) * bm_sapl_2D(ipts,ivm,1,iroot,icarbon) + fcn_wood(ivm) * & |
|---|
| 1072 | (bm_sapl_2D(ipts,ivm,1,isapabove,icarbon) + & |
|---|
| 1073 | bm_sapl_2D(ipts,ivm,1,isapbelow,icarbon) + & |
|---|
| 1074 | bm_sapl_2D(ipts,ivm,1,icarbres,icarbon))) |
|---|
| 1075 | |
|---|
| 1076 | ! Allocate the nitrogen |
|---|
| 1077 | cn_leaf=cn_leaf_init_2D(ipts,ivm) |
|---|
| 1078 | cn_wood=cn_leaf/fcn_wood(ivm) |
|---|
| 1079 | cn_root=cn_leaf/fcn_root(ivm) |
|---|
| 1080 | |
|---|
| 1081 | !+++CHECK+++ |
|---|
| 1082 | ! See comments at a similmar block of code for the |
|---|
| 1083 | ! tree saplings |
|---|
| 1084 | ! Use the C/N ratio to calculate the N content for |
|---|
| 1085 | ! all the biomass components. |
|---|
| 1086 | bm_sapl_2D(ipts,ivm,1,:,initrogen) = & |
|---|
| 1087 | bm_sapl_2D(ipts,ivm,1,:,icarbon) / cn_root |
|---|
| 1088 | |
|---|
| 1089 | ! Overwrite the N content for the leaves and wood |
|---|
| 1090 | ! components. Note that only sapabove is defined the |
|---|
| 1091 | ! other wood components are zero |
|---|
| 1092 | bm_sapl_2D(ipts,ivm,1,ileaf,initrogen) = & |
|---|
| 1093 | bm_sapl_2D(ipts,ivm,1,ileaf,icarbon) / cn_leaf |
|---|
| 1094 | bm_sapl_2D(ipts,ivm,1,isapabove,initrogen) = & |
|---|
| 1095 | bm_sapl_2D(ipts,ivm,1,isapabove,icarbon) / cn_wood |
|---|
| 1096 | !++++++++++++ |
|---|
| 1097 | |
|---|
| 1098 | ! Avoid deciduous PFTs to have leaves out at establishment |
|---|
| 1099 | ! Whether the saplings have leaves or don't have leaves the first |
|---|
| 1100 | ! year doesn't really matter. Either the approach is correct in the |
|---|
| 1101 | ! northern hemisphere or in the southern hemisphere. Anyhow, a |
|---|
| 1102 | ! spin-up is recommended to avoid issues with the initial conditions |
|---|
| 1103 | IF ( pheno_type(ivm) .NE. 1 ) THEN |
|---|
| 1104 | |
|---|
| 1105 | ! Not evergreen. Deciduous PFTs now need to survive half a year |
|---|
| 1106 | ! before bud burst. To ensure survival there are several options: |
|---|
| 1107 | ! (a) either the height of the initial vegetation is increased |
|---|
| 1108 | ! (this results in more reserves) or (b) the reserves could be |
|---|
| 1109 | ! increased. The second option may result in result in numerical |
|---|
| 1110 | ! issues further down the code as the optimal reserve level is |
|---|
| 1111 | ! calculated from the other biomass pools |
|---|
| 1112 | |
|---|
| 1113 | ! Notice that the sapwood is not put into reserves. Sapwood is |
|---|
| 1114 | ! important in phenology when the amount of leaves and roots |
|---|
| 1115 | ! are determined, so it needs to be already present. |
|---|
| 1116 | bm_sapl_2D(ipts,ivm,1,icarbres,:) = bm_sapl_2D(ipts,ivm,1,icarbres,:) + & |
|---|
| 1117 | bm_sapl_2D(ipts,ivm,1,ileaf,:) + & |
|---|
| 1118 | bm_sapl_2D(ipts,ivm,1,iroot,:) |
|---|
| 1119 | bm_sapl_2D(ipts,ivm,1,ileaf,:) = zero |
|---|
| 1120 | bm_sapl_2D(ipts,ivm,1,iroot,:) = zero |
|---|
| 1121 | |
|---|
| 1122 | ! When there are no leaves, the crop/grass is dormant |
|---|
| 1123 | plant_status(ipts,ivm) = idormant |
|---|
| 1124 | |
|---|
| 1125 | ELSE |
|---|
| 1126 | |
|---|
| 1127 | ! Initilize carbohydrate reserves for evergreen PFTs |
|---|
| 1128 | bm_sapl_2D(ipts,ivm,1,icarbres,:) = zero |
|---|
| 1129 | |
|---|
| 1130 | ! Evergreen plants always have a canopy |
|---|
| 1131 | plant_status(ipts,ivm) = icanopy |
|---|
| 1132 | |
|---|
| 1133 | ENDIF |
|---|
| 1134 | |
|---|
| 1135 | ! Debug |
|---|
| 1136 | IF (printlev_loc>=4) THEN |
|---|
| 1137 | DO iele=1,nelements |
|---|
| 1138 | WRITE(numout,*) ' young grass biomass (gC):',ivm,ipts |
|---|
| 1139 | WRITE(numout,*) ' bm_sapl_2D(ipts,:,1,ileaf,:,) ',& |
|---|
| 1140 | bm_sapl_2D(ipts,ivm,1,ileaf,iele) |
|---|
| 1141 | WRITE(numout,*) ' bm_sapl_2D(ipts,:,1,ispabove,:) ',& |
|---|
| 1142 | bm_sapl_2D(ipts,ivm,1,isapabove,iele) |
|---|
| 1143 | WRITE(numout,*) ' bm_sapl_2D(ipts,:,1,isapbelow,:) ',& |
|---|
| 1144 | bm_sapl_2D(ipts,ivm,1,isapbelow,iele) |
|---|
| 1145 | WRITE(numout,*) ' bm_sapl_2D(ipts,:,1,iheartabove,:) ',& |
|---|
| 1146 | bm_sapl_2D(ipts,ivm,1,iheartabove,iele) |
|---|
| 1147 | WRITE(numout,*) ' bm_sapl_2D(ipts,:,1,iheartbelow,:) ',& |
|---|
| 1148 | bm_sapl_2D(ipts,ivm,1,iheartbelow,iele) |
|---|
| 1149 | WRITE(numout,*) ' bm_sapl_2D(ipts,:,1,iroot,:) ',& |
|---|
| 1150 | bm_sapl_2D(ipts,ivm,1,iroot,iele) |
|---|
| 1151 | WRITE(numout,*) ' bm_sapl_2D(ipts,:,1,ifruit,:)', & |
|---|
| 1152 | bm_sapl_2D(ipts,ivm,1,ifruit,iele) |
|---|
| 1153 | WRITE(numout,*) ' bm_sapl_2D(ipts,:,1,icarbres,:)', & |
|---|
| 1154 | bm_sapl_2D(ipts,ivm,1,icarbres,iele) |
|---|
| 1155 | WRITE(numout,*) ' bm_sapl_2D(ipts,:,1,ilabile,:)', & |
|---|
| 1156 | bm_sapl_2D(ipts,ivm,1,ilabile,iele) |
|---|
| 1157 | ENDDO |
|---|
| 1158 | ENDIF |
|---|
| 1159 | ! - |
|---|
| 1160 | |
|---|
| 1161 | ! circ_class_biomass (gC tree-1). This |
|---|
| 1162 | ! allows to more easily check the mass balance closure |
|---|
| 1163 | circ_class_biomass(ipts,ivm,1,:,:) = bm_sapl_2D(ipts,ivm,1,:,:) |
|---|
| 1164 | circ_class_n(ipts,ivm,1) = tmp_ind |
|---|
| 1165 | |
|---|
| 1166 | ! Set leaf age classes -> all leaves will be current year leaves |
|---|
| 1167 | leaf_frac(ipts,ivm,:) = zero |
|---|
| 1168 | leaf_frac(ipts,ivm,1) = un |
|---|
| 1169 | |
|---|
| 1170 | ! Set time since last beginning of growing season but only |
|---|
| 1171 | ! for the first day of the whole simulation. When the model |
|---|
| 1172 | ! is initialized when_growthinit is set to undef. In subsequent |
|---|
| 1173 | ! time steps it should have a value. |
|---|
| 1174 | IF (when_growthinit(ipts,ivm) .EQ. undef) THEN |
|---|
| 1175 | when_growthinit(ipts,ivm) = 200 |
|---|
| 1176 | ENDIF |
|---|
| 1177 | |
|---|
| 1178 | ! The carbon and nitrogen biomass to build the saplings is taken |
|---|
| 1179 | ! from the atmosphere, keep track of amount to calculate the mass |
|---|
| 1180 | ! balance closure |
|---|
| 1181 | DO iele = 1,nelements |
|---|
| 1182 | atm_to_bm(ipts,ivm,iele) = atm_to_bm(ipts,ivm,iele) + & |
|---|
| 1183 | ((SUM(circ_class_biomass(ipts,ivm,1,:,iele))*& |
|---|
| 1184 | circ_class_n(ipts,ivm,1)) / dt) |
|---|
| 1185 | ENDDO |
|---|
| 1186 | |
|---|
| 1187 | ! Debug |
|---|
| 1188 | IF (printlev_loc>=4) THEN |
|---|
| 1189 | DO iele=1,nelements |
|---|
| 1190 | WRITE(numout,*) ' root to sapwood tradeoff (LF) : ', c0_alloc(ivm) |
|---|
| 1191 | WRITE(numout,*) ' grass sapling biomass: ',bm_sapl_2D(ipts,ivm,1,:,iele) |
|---|
| 1192 | ENDDO |
|---|
| 1193 | ENDIF |
|---|
| 1194 | ! - |
|---|
| 1195 | |
|---|
| 1196 | ELSEIF (.NOT. natural(ivm)) THEN |
|---|
| 1197 | |
|---|
| 1198 | ! Initialize croplands - else leaves_begin will never become true |
|---|
| 1199 | ! Set time since last beginning of growing season but only |
|---|
| 1200 | ! for the first day of the whole simulation. When the model |
|---|
| 1201 | ! is initialized when_growthinit is set to undef. In subsequent |
|---|
| 1202 | ! time steps it should have a value. |
|---|
| 1203 | IF (when_growthinit(ipts,ivm) .EQ. undef) THEN |
|---|
| 1204 | when_growthinit(ipts,ivm) = 200 |
|---|
| 1205 | ENDIF |
|---|
| 1206 | |
|---|
| 1207 | ! Debug |
|---|
| 1208 | IF (printlev_loc>=4) WRITE(numout,*) 'Initialize crops',ivm |
|---|
| 1209 | !- |
|---|
| 1210 | |
|---|
| 1211 | ENDIF ! .NOT. tree(ivm) |
|---|
| 1212 | |
|---|
| 1213 | ENDIF ! establish (for .NOT. is_tree(ivm)) |
|---|
| 1214 | |
|---|
| 1215 | !! 2.3.3 Declare PFT present |
|---|
| 1216 | ! Now that the PFT has biomass it should be declared 'present' |
|---|
| 1217 | ! everywhere in that grid box. Assign some additional properties |
|---|
| 1218 | IF (establish) THEN |
|---|
| 1219 | PFTpresent(ipts,ivm) = .TRUE. |
|---|
| 1220 | everywhere(ipts,ivm) = un |
|---|
| 1221 | age(ipts,ivm) = zero |
|---|
| 1222 | npp_longterm(ipts,ivm) = npp_longterm_init |
|---|
| 1223 | lm_lastyearmax(ipts,ivm) = SUM(bm_sapl_2D(ipts,ivm,:,ileaf,icarbon) * & |
|---|
| 1224 | circ_class_n(ipts,ivm,:)) |
|---|
| 1225 | ENDIF |
|---|
| 1226 | |
|---|
| 1227 | ENDIF ! establish or recruite |
|---|
| 1228 | |
|---|
| 1229 | ENDDO ! loop over pixels |
|---|
| 1230 | |
|---|
| 1231 | ENDDO ! loop over PFTs |
|---|
| 1232 | |
|---|
| 1233 | !! 3. Check numerical consistency of this routine |
|---|
| 1234 | |
|---|
| 1235 | IF (err_act.GT.1) THEN |
|---|
| 1236 | |
|---|
| 1237 | ! 3.1 Check surface area |
|---|
| 1238 | CALL check_vegetation_area("stomate_prescribe", npts, veget_max_begin, & |
|---|
| 1239 | veget_max,'pft') |
|---|
| 1240 | |
|---|
| 1241 | ! 3.2 Mass balance closure |
|---|
| 1242 | ! 3.2.1 Calculate final biomass |
|---|
| 1243 | pool_end(:,:,:) = zero |
|---|
| 1244 | DO ipar = 1,nparts |
|---|
| 1245 | DO iele = 1,nelements |
|---|
| 1246 | DO icir = 1,ncirc |
|---|
| 1247 | pool_end(:,:,iele) = pool_end(:,:,iele) + & |
|---|
| 1248 | (circ_class_biomass(:,:,icir,ipar,iele) * & |
|---|
| 1249 | circ_class_n(:,:,icir) * veget_max(:,:)) |
|---|
| 1250 | ENDDO |
|---|
| 1251 | ENDDO |
|---|
| 1252 | ENDDO |
|---|
| 1253 | |
|---|
| 1254 | ! 3.2.2 Calculate mass balance |
|---|
| 1255 | ! Common processes |
|---|
| 1256 | DO iele=1,nelements |
|---|
| 1257 | check_intern(:,:,iatm2land,iele) = check_intern(:,:,iatm2land,iele) + & |
|---|
| 1258 | atm_to_bm(:,:,iele) * veget_max(:,:) * dt |
|---|
| 1259 | check_intern(:,:,ipoolchange,iele) = -un * (pool_end(:,:,iele) - & |
|---|
| 1260 | pool_start(:,:,iele)) |
|---|
| 1261 | ENDDO |
|---|
| 1262 | |
|---|
| 1263 | closure_intern(:,:,:) = zero |
|---|
| 1264 | DO imbc = 1,nmbcomp |
|---|
| 1265 | DO iele=1,nelements |
|---|
| 1266 | ! Debug |
|---|
| 1267 | IF (printlev_loc>=3) THEN |
|---|
| 1268 | IF(iele .EQ. 1) WRITE(numout,"(9999(G12.5,:,','))") & |
|---|
| 1269 | 'imbc, test_pft, check_intern', test_grid , imbc, & |
|---|
| 1270 | test_pft, check_intern(test_grid,test_pft,imbc,iele) |
|---|
| 1271 | ENDIF |
|---|
| 1272 | !- |
|---|
| 1273 | closure_intern(:,:,iele) = closure_intern(:,:,iele) + & |
|---|
| 1274 | check_intern(:,:,imbc,iele) |
|---|
| 1275 | ENDDO |
|---|
| 1276 | ENDDO |
|---|
| 1277 | |
|---|
| 1278 | ! 3.3 Check mass balance closure |
|---|
| 1279 | CALL check_mass_balance("stomate_prescribe", closure_intern, npts, & |
|---|
| 1280 | pool_end, pool_start, veget_max, 'pft') |
|---|
| 1281 | |
|---|
| 1282 | ENDIF ! err_act .GT. 1 |
|---|
| 1283 | |
|---|
| 1284 | IF(printlev_loc>=3) WRITE(numout,*) 'Leaving stomate_prescribe.f90' |
|---|
| 1285 | |
|---|
| 1286 | END SUBROUTINE prescribe |
|---|
| 1287 | |
|---|
| 1288 | |
|---|
| 1289 | |
|---|
| 1290 | !! ================================================================================================================================ |
|---|
| 1291 | !! SUBROUTINE : make_sapling |
|---|
| 1292 | !! |
|---|
| 1293 | !>\BRIEF Given the dimensions and the allocation factors, calculate the biomass of a sapling |
|---|
| 1294 | !! |
|---|
| 1295 | !! DESCRIPTION (functional, design, flags): |
|---|
| 1296 | !! |
|---|
| 1297 | !! RECENT CHANGE(S): None |
|---|
| 1298 | !! |
|---|
| 1299 | !! MAIN OUTPUT VARIABLES(S): ::bm_sapl, ::KF |
|---|
| 1300 | !! |
|---|
| 1301 | !! REFERENCES : |
|---|
| 1302 | !! |
|---|
| 1303 | !! FLOWCHART : None |
|---|
| 1304 | !! \n |
|---|
| 1305 | !_ ================================================================================================================================ |
|---|
| 1306 | |
|---|
| 1307 | SUBROUTINE make_saplings(bm_sapl_2D_loc, tree_dia, & |
|---|
| 1308 | icir, ivm, ipts, KF_loc, c0_alloc_loc, plant_status_loc, & |
|---|
| 1309 | establish, cn_leaf_init_2D_loc,bypass_lm_loc) |
|---|
| 1310 | |
|---|
| 1311 | !! 0. Parameters and variables declaration |
|---|
| 1312 | |
|---|
| 1313 | !! 0.1 Input variables |
|---|
| 1314 | REAL(r_std), INTENT(in) :: tree_dia !! Tree diameter (m) |
|---|
| 1315 | INTEGER(i_std), INTENT(in) :: ivm, icir, ipts !! index (unitless) |
|---|
| 1316 | REAL(r_std), INTENT(in) :: c0_alloc_loc !! Root to sapwood tradeoff parameter |
|---|
| 1317 | LOGICAL, INTENT(in) :: establish !! yes: establish a new young vegetation. |
|---|
| 1318 | !! No: add saplings under an existing vegetation |
|---|
| 1319 | REAL(r_std), INTENT(in) :: cn_leaf_init_2D_loc!! initial leaf C/N ratio |
|---|
| 1320 | REAL(r_std), INTENT(in) :: KF_loc !! Scaling factor to convert sapwood mass |
|---|
| 1321 | !! into leaf mass (m) - this variable is |
|---|
| 1322 | !! passed to other routines |
|---|
| 1323 | |
|---|
| 1324 | !! 0.2 Output variables |
|---|
| 1325 | REAL(r_std), INTENT(out) :: bypass_lm_loc !! leaf mass before it is set to zero |
|---|
| 1326 | !! for deciduous PFTs. This value is |
|---|
| 1327 | !! in the calculation of KF (gC tree-1) |
|---|
| 1328 | !! 0.3 Modified variables |
|---|
| 1329 | REAL(r_std), DIMENSION(:,:),INTENT(inout) :: bm_sapl_2D_loc !! Sapling biomass for the functional |
|---|
| 1330 | REAL, INTENT(inout) :: plant_status_loc !! Growth and phenological status of the plant |
|---|
| 1331 | !! Different stati defined in constants |
|---|
| 1332 | !! 0.4 Local variables |
|---|
| 1333 | INTEGER(i_std) :: iele !! index (unitless) |
|---|
| 1334 | REAL(r_std) :: cn_leaf,cn_root !! CN ratio of leaves and roots (gC/gN) |
|---|
| 1335 | REAL(r_std) :: cn_wood !! CN ratio of wood pool (gC/gN) |
|---|
| 1336 | REAL(r_std) :: tree_height !! Tree height (m) |
|---|
| 1337 | |
|---|
| 1338 | !_ ================================================================================================================================ |
|---|
| 1339 | |
|---|
| 1340 | ! The woody biomass is contained in four components. Thus, |
|---|
| 1341 | ! wood_init = isapabove + isapbelow + iheartabove + iheartbelow. |
|---|
| 1342 | ! Given that isapbelow = isapbelow and iheartabove = |
|---|
| 1343 | ! iheartbelow = bm_sapl_heartabove*isapabove. If |
|---|
| 1344 | ! bm_sapl_heartbelow = bm_sapl_heartabove = 0.2, then |
|---|
| 1345 | ! isapabove = wood_init/2.4 |
|---|
| 1346 | |
|---|
| 1347 | IF ((dia_init_max(ivm)+dia_init_min(ivm))/2 .GT. 0.05) THEN |
|---|
| 1348 | CALL ipslerr_p(3,'Initial diameter of sapling was set to more than 5cm',& |
|---|
| 1349 | 'The current estimate of a sapling biomass does not account & |
|---|
| 1350 | for heartwood','If you want to use larger diameters add a & |
|---|
| 1351 | heartwood estimate in the code','') |
|---|
| 1352 | ! bm_sapl_2D_loc(isapabove,icarbon) = |
|---|
| 1353 | ! tree_dia**(2+pipe_tune3(ivm))*pi/4* & |
|---|
| 1354 | ! 2*tree_ff(ivm)*pipe_density(ivm)*pipe_tune2(ivm) / & |
|---|
| 1355 | ! (2. + bm_sapl_heartabove + bm_sapl_heartbelow) |
|---|
| 1356 | !bm_sapl_2D_loc(isapbelow,icarbon) = bm_sapl_2D_loc(isapabove,icarbon) |
|---|
| 1357 | !bm_sapl_2D_loc(iheartabove,icarbon) = bm_sapl_heartabove * & |
|---|
| 1358 | ! bm_sapl_2D_loc(isapabove,icarbon) |
|---|
| 1359 | !bm_sapl_2D_loc(iheartbelow,icarbon) = bm_sapl_heartbelow * & |
|---|
| 1360 | ! bm_sapl_2D_loc(isapbelow,icarbon) |
|---|
| 1361 | |
|---|
| 1362 | ELSE |
|---|
| 1363 | |
|---|
| 1364 | bm_sapl_2D_loc(iheartabove,icarbon) = 0.0001 |
|---|
| 1365 | bm_sapl_2D_loc(iheartbelow,icarbon) = 0.0001 |
|---|
| 1366 | bm_sapl_2D_loc(isapabove,icarbon) = tree_dia**(2+pipe_tune3(ivm))*pi/4* & |
|---|
| 1367 | tree_ff(ivm)*pipe_density(ivm)*pipe_tune2(ivm) - & |
|---|
| 1368 | bm_sapl_2D_loc(iheartabove,icarbon) |
|---|
| 1369 | bm_sapl_2D_loc(isapbelow,icarbon) = bm_sapl_2D_loc(isapabove,icarbon) |
|---|
| 1370 | |
|---|
| 1371 | ENDIF |
|---|
| 1372 | |
|---|
| 1373 | ! Use the allometric relationships to calculate initial leaf |
|---|
| 1374 | ! and root mass. Note that wstress should be accounted for |
|---|
| 1375 | ! through c0_alloc. There was an inconsistency in the |
|---|
| 1376 | ! original calculation (OCN) - most pools were in gN but |
|---|
| 1377 | ! leaves were in gC. A correction is proposed. Given that |
|---|
| 1378 | ! icarbes was just set to zero the equation for ilabile is |
|---|
| 1379 | ! a bit overkill and icarbes could be omitted. |
|---|
| 1380 | tree_height = ((((bm_sapl_2D_loc(isapabove,icarbon)+ & |
|---|
| 1381 | bm_sapl_2D_loc(isapbelow,icarbon)+bm_sapl_2D_loc(iheartabove,icarbon)+ & |
|---|
| 1382 | bm_sapl_2D_loc(iheartbelow,icarbon))/(tree_ff(ivm)* & |
|---|
| 1383 | pipe_density(ivm))*4/pi)**(pipe_tune3(ivm)/2)) & |
|---|
| 1384 | *pipe_tune2(ivm))**(1/(pipe_tune3(ivm)/2+1)) |
|---|
| 1385 | |
|---|
| 1386 | bm_sapl_2D_loc(ileaf,icarbon) = ( bm_sapl_2D_loc(isapabove,icarbon) + & |
|---|
| 1387 | bm_sapl_2D_loc(isapbelow,icarbon) ) * KF_loc / tree_height |
|---|
| 1388 | bm_sapl_2D_loc(iroot,icarbon) = bm_sapl_2D_loc(ileaf,icarbon) / & |
|---|
| 1389 | ( KF_loc * c0_alloc_loc ) |
|---|
| 1390 | bm_sapl_2D_loc(icarbres,icarbon)=zero |
|---|
| 1391 | bm_sapl_2D_loc(ifruit,icarbon) = zero |
|---|
| 1392 | bm_sapl_2D_loc(ilabile,icarbon) = labile_to_total * & |
|---|
| 1393 | (bm_sapl_2D_loc(ileaf,icarbon) + & |
|---|
| 1394 | fcn_root(ivm) * bm_sapl_2D_loc(iroot,icarbon) + & |
|---|
| 1395 | fcn_wood(ivm) * (bm_sapl_2D_loc(isapabove,icarbon) + & |
|---|
| 1396 | bm_sapl_2D_loc(isapbelow,icarbon) + & |
|---|
| 1397 | bm_sapl_2D_loc(icarbres,icarbon))) |
|---|
| 1398 | |
|---|
| 1399 | ! Debug |
|---|
| 1400 | IF (printlev_loc>=4) THEN |
|---|
| 1401 | WRITE(numout,*) ' +++++ CHECK INSIDE make_sapling ++++ ' |
|---|
| 1402 | WRITE(numout,*) 'Circumference class: ',icir,' PFT type: ',ivm |
|---|
| 1403 | WRITE(numout,*) ' root to sapwood tradeoff p :', c0_alloc_loc |
|---|
| 1404 | WRITE(numout,*) ' sapling height, sapling diameter ', & |
|---|
| 1405 | tree_height , ( tree_height/pipe_tune2(ivm) ) ** & |
|---|
| 1406 | ( 1. / pipe_tune3(ivm) ) |
|---|
| 1407 | WRITE(numout,*) 'pipe_density, ',pipe_density(ivm) |
|---|
| 1408 | ENDIF |
|---|
| 1409 | ! - |
|---|
| 1410 | |
|---|
| 1411 | ! Allocate the nitrogen |
|---|
| 1412 | cn_leaf=cn_leaf_init_2D_loc |
|---|
| 1413 | cn_wood=cn_leaf/fcn_wood(ivm) |
|---|
| 1414 | cn_root=cn_leaf/fcn_root(ivm) |
|---|
| 1415 | |
|---|
| 1416 | ! Use the C/N ratio to calculate the N content for |
|---|
| 1417 | ! all the biomass components. |
|---|
| 1418 | bm_sapl_2D_loc(:,initrogen) = & |
|---|
| 1419 | bm_sapl_2D_loc(:,icarbon) / cn_root |
|---|
| 1420 | |
|---|
| 1421 | ! Overwrite the N content for the leaves and wood |
|---|
| 1422 | ! components. |
|---|
| 1423 | bm_sapl_2D_loc(ileaf,initrogen) = & |
|---|
| 1424 | bm_sapl_2D_loc(ileaf,icarbon) / cn_leaf |
|---|
| 1425 | bm_sapl_2D_loc(isapabove,initrogen) = & |
|---|
| 1426 | bm_sapl_2D_loc(isapabove,icarbon) / cn_wood |
|---|
| 1427 | bm_sapl_2D_loc(isapbelow,initrogen) = & |
|---|
| 1428 | bm_sapl_2D_loc(isapbelow,icarbon) / cn_wood |
|---|
| 1429 | bm_sapl_2D_loc(iheartabove,initrogen) = & |
|---|
| 1430 | bm_sapl_2D_loc(iheartabove,icarbon) / cn_wood |
|---|
| 1431 | bm_sapl_2D_loc(iheartbelow,initrogen) = & |
|---|
| 1432 | bm_sapl_2D_loc(iheartbelow,icarbon) / cn_wood |
|---|
| 1433 | |
|---|
| 1434 | ! Store the leaf mass so it can be used in the calculation of KF |
|---|
| 1435 | ! when using a dynamic sla. |
|---|
| 1436 | bypass_lm_loc = bm_sapl_2D_loc(ileaf,icarbon) |
|---|
| 1437 | |
|---|
| 1438 | ! Avoid deciduous PFTs to have leaves out at establishment |
|---|
| 1439 | ! Whether the saplings have leaves or don't have leaves the |
|---|
| 1440 | ! first year doesn't really matter. Either the approach is |
|---|
| 1441 | ! correct in the northern hemisphere or in the southern |
|---|
| 1442 | ! hemisphere. Anyhow, a spin-up is needed to avoid issues |
|---|
| 1443 | ! with the initial conditions |
|---|
| 1444 | IF ( pheno_type(ivm) .NE. 1 ) THEN |
|---|
| 1445 | |
|---|
| 1446 | ! Not evergreen. Deciduous PFTs now need to survive an extra |
|---|
| 1447 | ! year before bud burst. To ensure survival there are several |
|---|
| 1448 | ! options: (a) either the height of the initial vegetation is |
|---|
| 1449 | ! increased (this results in more reserves) or (b) the reserves |
|---|
| 1450 | ! could be increased. The second option may result in numerical |
|---|
| 1451 | ! issues further down the code as the optimal reserve level is |
|---|
| 1452 | ! calculated from the other biomass pools. Also some initial |
|---|
| 1453 | ! tests showed that higher reserves simply resulted in more |
|---|
| 1454 | ! respiration. |
|---|
| 1455 | ! Tip: use taller trees to start with if they die between |
|---|
| 1456 | ! prescribe and leaf onset |
|---|
| 1457 | bm_sapl_2D_loc(icarbres,:) = (bm_sapl_2D_loc(icarbres,:) + & |
|---|
| 1458 | bm_sapl_2D_loc(ileaf,:) + bm_sapl_2D_loc(iroot,:)) |
|---|
| 1459 | bm_sapl_2D_loc(ileaf,:) = zero |
|---|
| 1460 | bm_sapl_2D_loc(iroot,:) = zero |
|---|
| 1461 | |
|---|
| 1462 | ! When deciduous trees have no leaves they are dormant |
|---|
| 1463 | ! If the code is used for recruitment use the actual |
|---|
| 1464 | ! plant_status of the overstorey trees for the saplings |
|---|
| 1465 | IF (establish) THEN |
|---|
| 1466 | plant_status_loc = idormant |
|---|
| 1467 | ENDIF |
|---|
| 1468 | |
|---|
| 1469 | ELSE |
|---|
| 1470 | |
|---|
| 1471 | ! Initilize carbohydrate reserves for evergreen PFTs |
|---|
| 1472 | bm_sapl_2D_loc(icarbres,:) = zero |
|---|
| 1473 | |
|---|
| 1474 | ! Evergreen trees always have a canopy. If the code |
|---|
| 1475 | ! is used for recruitment use the actual plant_status |
|---|
| 1476 | ! of the overstorey trees for the saplings |
|---|
| 1477 | IF (establish) THEN |
|---|
| 1478 | plant_status_loc = icanopy |
|---|
| 1479 | ENDIF |
|---|
| 1480 | |
|---|
| 1481 | ENDIF |
|---|
| 1482 | |
|---|
| 1483 | ! Debug |
|---|
| 1484 | IF (printlev_loc>=4) THEN |
|---|
| 1485 | DO iele=1,nelements |
|---|
| 1486 | ! This is the only use of icir, ivm, and ipts in this routine. |
|---|
| 1487 | ! This could be placed outside the routine, but it's a bit |
|---|
| 1488 | ! cleaner here. |
|---|
| 1489 | WRITE(numout,*) ' sapling biomass (gC):',icir,ivm,ipts |
|---|
| 1490 | WRITE(numout,*) ' bm_sapl_2D_loc(ileaf,:,) ',& |
|---|
| 1491 | bm_sapl_2D_loc(ileaf,iele) |
|---|
| 1492 | WRITE(numout,*) ' bm_sapl_2D_loc(ispabove,:) ',& |
|---|
| 1493 | bm_sapl_2D_loc(isapabove,iele) |
|---|
| 1494 | WRITE(numout,*) ' bm_sapl_2D_loc(isapbelow,:) ',& |
|---|
| 1495 | bm_sapl_2D_loc(isapbelow,iele) |
|---|
| 1496 | WRITE(numout,*) ' bm_sapl_2D_loc(iheartabove,:) ',& |
|---|
| 1497 | bm_sapl_2D_loc(iheartabove,iele) |
|---|
| 1498 | WRITE(numout,*) ' bm_sapl_2D_loc(iheartbelow,:) ',& |
|---|
| 1499 | bm_sapl_2D_loc(iheartbelow,iele) |
|---|
| 1500 | WRITE(numout,*) ' bm_sapl_2D_loc(iroot,:) ',& |
|---|
| 1501 | bm_sapl_2D_loc(iroot,iele) |
|---|
| 1502 | WRITE(numout,*) ' bm_sapl_2D_loc(ifruit,:)', & |
|---|
| 1503 | bm_sapl_2D_loc(ifruit,iele) |
|---|
| 1504 | WRITE(numout,*) ' bm_sapl_2D_loc(icarbres,:)', & |
|---|
| 1505 | bm_sapl_2D_loc(icarbres,iele) |
|---|
| 1506 | WRITE(numout,*) ' bm_sapl_2D_loc(ilabile,:)', & |
|---|
| 1507 | bm_sapl_2D_loc(ilabile,iele) |
|---|
| 1508 | ENDDO |
|---|
| 1509 | ENDIF |
|---|
| 1510 | |
|---|
| 1511 | END SUBROUTINE make_saplings |
|---|
| 1512 | |
|---|
| 1513 | END MODULE stomate_prescribe |
|---|
