mode_cotwo.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     ##################
      MODULE mode_cotwo
!     ##################
!
!!****  *MODE_COTWO * - contains some needed computations for vegetation
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!    
!!    EXTERNAL
!!    --------
!!       
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!      A. Boone                 * Meteo France *       
!!
!!    MODIFICATIONS
!!    -------------
!!      Original        05/03/15
!-----------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!
!
INTERFACE gauleg
   MODULE PROCEDURE gauleg
END INTERFACE
!
!
!-------------------------------------------------------------------------------
 CONTAINS
!
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE gauleg(PX1,PX2,PX,PW,KN)
!
!
!!****  *GAULEG*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the Gaussian weights for integration of net assimilation
!     and stomatal conductance over the canopy depth
!         
!     
!!**  METHOD
!!    ------
!
!     1) Calculate the weights and abscissa using Gaussian Quadrature
!
!!    EXTERNAL
!!    --------
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!    MODD_CST
!!      
!!    REFERENCE
!!    ---------
!!
!!    Calvet et al. (1998) For. Agri. Met.
!!      
!!    AUTHOR
!!    ------
!!
!!      A. Boone           * Meteo-France *
!!      (following Belair)
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    27/10/97 
!!
!-------------------------------------------------------------------------------
!
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_csts, ONLY : xpi
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
INTEGER,             INTENT(IN)   :: kn
!                                    number of points at which Gaussian
!                                    weights are evaluated/needed
!
REAL,                INTENT(IN)   :: px1, px2
!                                    mathematical/numerical values needed for 
!                                    weight computation
!
REAL, DIMENSION(KN), INTENT(OUT)  :: px, pw
!                                    PX = abscissa
!                                    PW = Gaussian weights
!
!
!*      0.2    declarations of local variables
!
REAL, PARAMETER                             :: ppeps = 3.0e-14
!                                              convergence criteria
!
INTEGER ji,jj,jk                             ! loop indexes
!
INTEGER im                                   ! 
!
REAL zxm, zxl, zz, zp1, zp2, zp3, zpp, zz1   ! dummy variables needed for 
!
REAL(KIND=JPRB) :: zhook_handle
!                                              computation of the gaussian weights
!
!-------------------------------------------------------------------------------
!
!  calcul des poids et abscisses par la methode de quad de Gauss
!
IF (lhook) CALL dr_hook('GAULEG',0,zhook_handle)
!
im  = (kn+1)/2
zxm = 0.50*(px2+px1)
zxl = 0.50*(px2-px1)
!      
im_point_loop: DO ji = 1,im
   zz  = cos(xpi*(float(ji)-.250)/(float(kn)+.50))
!
!  begin iteration:
!
   iteration_loop: DO jk = 1,100
      zp1 = 1.
      zp2 = 0.
      DO jj = 1,kn
         zp3 = zp2
         zp2 = zp1
         zp1 = ((2.*(jj)-1.)*zz*zp2-(float(jj)-1.)*zp3)/jj
      END DO
      zpp = float(kn)*(zz*zp1-zp2)/(zz*zz-1.)
      zz1 = zz
      zz  = zz1-zp1/zpp
      IF(abs(zz-zz1).LE.ppeps)EXIT
   END DO iteration_loop
!
!  end iteration.
!
!  compute abscissa
!
   px(ji)      = zxm - zxl*zz
   px(kn+1-ji) = zxm + zxl*zz
!
!  compute weights
!
   pw(ji)      = 2.0*zxl/((1.0-zz*zz)*zpp*zpp)
   pw(kn+1-ji) = pw(ji)
END DO im_point_loop
!
IF (lhook) CALL dr_hook('GAULEG',1,zhook_handle)
!
!
END SUBROUTINE gauleg
!####################################################################
!####################################################################
!####################################################################
!
END MODULE mode_cotwo