source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_stomate/stomate_vmax.f90 @ 880

! calculates the leaf efficiency
!
!< $HeadURL$ 
!< $Date$
!< $Author$
!< $Revision$
! IPSL (2006)
!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
MODULE stomate_vmax

  ! modules used:

  USE ioipsl
  USE stomate_data
  USE constantes
  USE pft_parameters

  IMPLICIT NONE

  ! private & public routines

  PRIVATE
  PUBLIC vmax, vmax_clear

  ! first call
  LOGICAL, SAVE                                              :: firstcall = .TRUE.

CONTAINS

  SUBROUTINE vmax_clear
    firstcall=.TRUE.
  END SUBROUTINE vmax_clear

  SUBROUTINE vmax (npts, dt, &
       leaf_age, leaf_frac, &
       vcmax, vjmax)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                        :: npts
    ! time step of Stomate in days
    REAL(r_std), INTENT(in)                                     :: dt

    ! 0.2 modified fields

    ! leaf age (days)
    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_age
    ! fraction of leaves in leaf age class
    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_frac

    ! 0.3 output

    ! Maximum rate of carboxylation
    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: vcmax
    ! Maximum rate of RUbp regeneration
    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: vjmax

    ! 0.4 local

    ! leaf efficiency (vcmax/vcmax_opt)
    REAL(r_std), DIMENSION(npts)                                :: leaf_efficiency
    ! change of fraction of leaves in age class
    REAL(r_std), DIMENSION(npts,nvm,nleafages)                 :: d_leaf_frac
    ! new leaf age (d)
    REAL(r_std), DIMENSION(npts,nleafages)                      :: leaf_age_new
    ! sum of leaf age fractions, for normalization
    REAL(r_std), DIMENSION(npts)                                :: sumfrac
    ! relative leaf age (age/critical age)
    REAL(r_std), DIMENSION(npts)                                :: rel_age
    ! Index
    INTEGER(i_std)                                              :: j,m

    ! =========================================================================

    IF (bavard.GE.3) WRITE(numout,*) 'Entering vmax'

    !
    ! 1 Initialization
    !

    !
    ! 1.1 first call: info about flags and parameters.
    !

    IF ( firstcall ) THEN

       WRITE(numout,*) 'vmax:'

       WRITE(numout,*) '   > offset (minimum vcmax/vmax_opt):' , vmax_offset
       WRITE(numout,*) '   > relative leaf age at which vmax attains vcmax_opt:', leafage_firstmax
       WRITE(numout,*) '   > relative leaf age at which vmax falls below vcmax_opt:', leafage_lastmax
       WRITE(numout,*) '   > relative leaf age at which vmax attains its minimum:', leafage_old

       firstcall = .FALSE.

    ENDIF

    !
    ! 1.2 initialize output
    !

    vcmax(:,:) = zero
    vjmax(:,:) = zero

    !
    ! 2 leaf age: general increase and turnover between age classes.
    !

    !
    ! 2.1 increase leaf age
    !

    DO m = 1, nleafages

       DO j = 2,nvm
          WHERE ( leaf_frac(:,j,m) .GT. min_stomate )

             leaf_age(:,j,m) = leaf_age(:,j,m) + dt
             
          ENDWHERE
       ENDDO

    ENDDO

    !
    ! 2.2 turnover between leaf age classes
    !     d_leaf_frac(:,:,m) = what leaves m-1 and goes into m
    !

    DO j = 2,nvm

       ! 2.2.1 fluxes

       ! nothing goes into first age class
       d_leaf_frac(:,j,1) = zero

       ! from m-1 to m
       DO m = 2, nleafages 

          d_leaf_frac(:,j,m) = leaf_frac(:,j,m-1) * dt/leaf_timecst(j)

       ENDDO

       ! 2.2.2 new leaf age in class
       !       new age = ( old age * old fraction + fractional increase * age of source ) /
       !                 new fraction

       leaf_age_new(:,:) = zero

       DO m = 2, nleafages-1
          !      DO m=2, nleafages

          WHERE ( d_leaf_frac(:,j,m) .GT. min_stomate )

             leaf_age_new(:,m) = ( ( (leaf_frac(:,j,m)- d_leaf_frac(:,j,m+1)) * leaf_age(:,j,m) )  + &
                  ( d_leaf_frac(:,j,m) * leaf_age(:,j,m-1) ) ) / &
                  ( leaf_frac(:,j,m) + d_leaf_frac(:,j,m)- d_leaf_frac(:,j,m+1) )

             !           leaf_age_new(:,m) = ( ( leaf_frac(:,j,m) * leaf_age(:,j,m) )  + &
             !                                ( d_leaf_frac(:,j,m) * leaf_age(:,j,m-1) ) ) / &
             !                              ( leaf_frac(:,j,m) + d_leaf_frac(:,j,m) )

          ENDWHERE

       ENDDO       ! Loop over age classes

       WHERE ( d_leaf_frac(:,j,nleafages) .GT. min_stomate )

          leaf_age_new(:,nleafages) = ( ( leaf_frac(:,j,nleafages) * leaf_age(:,j,nleafages) )  + &
               ( d_leaf_frac(:,j,nleafages) * leaf_age(:,j,nleafages-1) ) ) / &
               ( leaf_frac(:,j,nleafages) + d_leaf_frac(:,j,nleafages) )

       ENDWHERE

       DO m = 2, nleafages

          WHERE ( d_leaf_frac(:,j,m) .GT. min_stomate )

             leaf_age(:,j,m) = leaf_age_new(:,m)

          ENDWHERE

       ENDDO       ! Loop over age classes

       ! 2.2.3 calculate new fraction

       DO m = 2, nleafages

          ! where the change comes from
          leaf_frac(:,j,m-1) = leaf_frac(:,j,m-1) - d_leaf_frac(:,j,m)

          ! where it goes to
          leaf_frac(:,j,m) = leaf_frac(:,j,m) + d_leaf_frac(:,j,m)

       ENDDO

       ! 2.2.4 renormalize fractions in order to prevent accumulation 
       !       of numerical errors

       ! correct small negative values

       DO m = 1, nleafages
          leaf_frac(:,j,m) = MAX( zero, leaf_frac(:,j,m) )
       ENDDO

       ! total of fractions, should be very close to one where there is leaf mass

       sumfrac(:) = zero

       DO m = 1, nleafages

          sumfrac(:) = sumfrac(:) + leaf_frac(:,j,m)

       ENDDO

       ! normalize

       DO m = 1, nleafages

          WHERE ( sumfrac(:) .GT. min_stomate )

             leaf_frac(:,j,m) = leaf_frac(:,j,m) / sumfrac(:) 

          ELSEWHERE

             leaf_frac(:,j,m) = zero

          ENDWHERE

       ENDDO

    ENDDO         ! Loop over PFTs

    !
    ! 3 calculate vmax as a function of the age
    !

    DO j = 2,nvm

       vcmax(:,j) = zero
       vjmax(:,j) = zero

       ! sum up over the different age classes

       DO m = 1, nleafages

          !
          ! 3.1 efficiency in each of the age classes
          !     increases from 0 to 1 at the beginning (rel_age < leafage_firstmax), stays 1
          !     until rel_age = leafage_lastmax, then decreases to vmax_offset at 
          !     rel_age = leafage_old, then stays at vmax_offset.
          !

          rel_age(:) = leaf_age(:,j,m) / leafagecrit(j)

          leaf_efficiency(:) = MAX( vmax_offset, MIN( un, &
               vmax_offset + (un - vmax_offset) * rel_age(:) / leafage_firstmax, &
               un - (un - vmax_offset) * ( rel_age(:) - leafage_lastmax ) / &
               ( leafage_old - leafage_lastmax ) ) )

          !
          ! 3.2 add to mean vmax
          !

          vcmax(:,j) = vcmax(:,j) + vcmax_opt(j) * leaf_efficiency(:) * leaf_frac(:,j,m)
          vjmax(:,j) = vjmax(:,j) + vjmax_opt(j) * leaf_efficiency(:) * leaf_frac(:,j,m)

       ENDDO     ! loop over age classes

    ENDDO       ! loop over PFTs

    IF (bavard.GE.4) WRITE(numout,*) 'Leaving vmax'

  END SUBROUTINE vmax

END MODULE stomate_vmax