av_pgd_param.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.
!     ################################################################
      SUBROUTINE av_pgd_param (DTI, &
                               pfield,pvegtype,pdata,hsftype,hatype,pdz,kdecade)
!     ################################################################
!
!!**** *AV_PATCH_PGD* average for each surface patch a secondary physiographic 
!!                    variable from the
!!              fractions of coverage class.
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!
!!    The averaging is performed with one way into three:
!!
!!    - arithmetic averaging (HATYPE='ARI')
!!
!!    - inverse    averaging (HATYPE='INV')
!!
!!    - inverse of square logarithm averaging (HATYPE='CDN') :
!!
!!      1 / ( ln (dz/data) )**2
!!
!!      This latest uses (if available) the height of the first model mass
!!      level. In the other case, 20m is chosen. It works for roughness lengths.
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    F.Solmon /V. Masson       
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    15/12/97
!!    V. Masson   01/2004  Externalization
!!    R. Alkama   04/2012  add 6 new tree vegtype (9 instead 3)
!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
!
USE modd_data_isba_n, ONLY : data_isba_t
!
USE modd_surf_par,       ONLY : xundef
USE modd_data_cover_par, ONLY : nvt_tebd, nvt_bone, nvt_trbe, nvt_trbd, nvt_tebe,  &
                                nvt_tene, nvt_bobd, nvt_bond, nvt_shrb, nvegtype,  &
                                xcdref

!
USE modi_vegtype_to_patch 
!
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
!
TYPE(data_isba_t), INTENT(INOUT) :: dti
!
REAL, DIMENSION(:,:), INTENT(OUT) :: pfield  ! secondary field to construct
REAL, DIMENSION(:,:), INTENT(IN)  :: pvegtype  ! fraction of each cover class
REAL, DIMENSION(:,:), INTENT(IN)  :: pdata   ! secondary field value for each class
 CHARACTER(LEN=3),     INTENT(IN)  :: hsftype ! Type of surface where the field
                                               ! is defined
 CHARACTER(LEN=3),     INTENT(IN)  :: hatype  ! Type of averaging
REAL, DIMENSION(:),   INTENT(IN), OPTIONAL :: pdz    ! first model half level
INTEGER,              INTENT(IN), OPTIONAL :: kdecade ! current month
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
!
INTEGER :: icover  ! number of cover classes
INTEGER :: jcover  ! loop on cover classes
!
! nbe of vegtype
! nbre of patches
INTEGER :: jvegtype! loop on vegtype
INTEGER :: ipatch  ! number of patches
INTEGER :: jpatch  ! PATCH index
INTEGER :: jj, ji
!
REAL, DIMENSION(SIZE(PFIELD,1),NVEGTYPE)  :: zweight
!
REAL, DIMENSION(SIZE(PFIELD,1),SIZE(PFIELD,2))   :: zsum_weight_patch
!
REAL, DIMENSION(SIZE(PFIELD,1),SIZE(PFIELD,2))   :: zwork
REAL, DIMENSION(SIZE(PFIELD,1),SIZE(PFIELD,2))   :: zdz
!
INTEGER, DIMENSION(SIZE(PFIELD,1),SIZE(PFIELD,2))  :: nmask
INTEGER, DIMENSION(SIZE(PFIELD,2)) :: jcount
INTEGER ::  patch_list(nvegtype)
REAL(KIND=JPRB) :: zhook_handle

!-------------------------------------------------------------------------------
!
!*    1.1    field does not exist
!            --------------------
!
IF (lhook) CALL dr_hook('AV_PGD_PARAM',0,zhook_handle)
IF (SIZE(pfield)==0 .AND. lhook) CALL dr_hook('AV_PGD_PARAM',1,zhook_handle)
IF (SIZE(pfield)==0) RETURN
!
!-------------------------------------------------------------------------------
!
!*    1.2    Initializations
!            ---------------
!
ipatch=SIZE(pfield,2)
!
!
IF (present(pdz)) THEN
  DO jpatch=1,ipatch
      zdz(:,jpatch)=pdz(:)
  END DO
ELSE
  zdz(:,:)=xcdref
END IF
!
pfield(:,:)=xundef
!
zwork(:,:)=0.
zweight(:,:)=0.
zsum_weight_patch(:,:)=0.
!
DO jvegtype=1,nvegtype
  patch_list(jvegtype) = vegtype_to_patch(jvegtype, ipatch)
ENDDO

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!*    2.     Selection of the weighting function for vegtype
!            -----------------------------------
!
SELECT CASE (hsftype)

   CASE('NAT','GRD')
     DO jvegtype=1,nvegtype
       zweight(:,jvegtype)=pvegtype(:,jvegtype)
     END DO

   CASE('VEG','GRV')
     DO jvegtype=1,nvegtype
       zweight(:,jvegtype)=pvegtype(:,jvegtype)*dti%XPAR_VEG(:,kdecade,jvegtype)
     END DO

   CASE('BAR','GRB')
     DO jvegtype=1,nvegtype
       zweight(:,jvegtype)=pvegtype(:,jvegtype)*(1.-dti%XPAR_VEG(:,kdecade,jvegtype))
     END DO
     
    CASE('DVG','GDV') ! for diffusion scheme only, average only on vegetated area
     DO jvegtype=1,nvegtype
       WHERE ( sum(dti%XPAR_LAI(:,:,jvegtype),2) .GT. 0.0) &
        zweight(:,jvegtype)=pvegtype(:,jvegtype)
     END DO

   CASE('LAI','GRL')
     DO jvegtype=4,nvegtype
       zweight(:,jvegtype)=pvegtype(:,jvegtype)*dti%XPAR_LAI(:,kdecade,jvegtype)
     END DO

    CASE('TRE','GRT')
      zweight(:,:)=0.
      WHERE (pvegtype(:,nvt_tebd)>0.)
        zweight(:,nvt_tebd)=pvegtype(:,nvt_tebd)
      ENDWHERE
      WHERE (pvegtype(:,nvt_bone)>0.)
        zweight(:,nvt_bone)=pvegtype(:,nvt_bone)
      ENDWHERE
      WHERE (pvegtype(:,nvt_trbe)>0.)
        zweight(:,nvt_trbe)=pvegtype(:,nvt_trbe)
      ENDWHERE

      WHERE (pvegtype(:,nvt_trbd)>0.)
        zweight(:,nvt_trbd)=pvegtype(:,nvt_trbd)
      ENDWHERE
      WHERE (pvegtype(:,nvt_tebe)>0.)
        zweight(:,nvt_tebe)=pvegtype(:,nvt_tebe)
      ENDWHERE
      WHERE (pvegtype(:,nvt_tene)>0.)
        zweight(:,nvt_tene)=pvegtype(:,nvt_tene)
      ENDWHERE
      WHERE (pvegtype(:,nvt_bobd)>0.)
        zweight(:,nvt_bobd)=pvegtype(:,nvt_bobd)
      ENDWHERE
      WHERE (pvegtype(:,nvt_bond)>0.)
        zweight(:,nvt_bond)=pvegtype(:,nvt_bond)
      ENDWHERE
      WHERE (pvegtype(:,nvt_shrb)>0.)
        zweight(:,nvt_shrb)=pvegtype(:,nvt_shrb)
      ENDWHERE

    CASE default
       CALL abor1_sfx('AV_PGD_PARAM: WEIGHTING FUNCTION FOR VEGTYPE NOT ALLOWED')
END SELECT
!
!-------------------------------------------------------------------------------
!
!*    3.     Averaging
!            ---------
!
!*    3.1    Work arrays given for each patch
!            -----------
!
!*    3.2    Selection of averaging type
!            ---------------------------
!
SELECT CASE (hatype)
!
!-------------------------------------------------------------------------------
!
!*    3.3    Arithmetic averaging
!            --------------------
!
  CASE ('ARI')
!
    DO jvegtype=1,nvegtype
      jpatch= patch_list(jvegtype)
      DO jj=1,SIZE(pdata,1)
        zsum_weight_patch(jj,jpatch) = zsum_weight_patch(jj,jpatch) + zweight(jj,jvegtype)
        zwork(jj,jpatch) =  zwork(jj,jpatch) + pdata(jj,jvegtype)  * zweight(jj,jvegtype)
      ENDDO
    END DO
!
!-------------------------------------------------------------------------------
!
!*    3.4    Inverse averaging
!            -----------------
!
  CASE('INV' )
!
   DO jvegtype=1,nvegtype 
     jpatch=patch_list(jvegtype) 
     DO jj=1,SIZE(pdata,1)
       zsum_weight_patch(jj,jpatch) = zsum_weight_patch(jj,jpatch)+zweight(jj,jvegtype)
       IF (pdata(jj,jvegtype).NE.0.) THEN
         zwork(jj,jpatch)= zwork(jj,jpatch) + 1./ pdata(jj,jvegtype) * zweight(jj,jvegtype)
       ENDIF
     ENDDO
   END DO
!
!-------------------------------------------------------------------------------!
!
!*    3.5    Roughness length averaging
!            --------------------------

!
  CASE('CDN')
!
    DO jvegtype=1,nvegtype
      jpatch=patch_list(jvegtype)
      DO jj=1,SIZE(pdata,1)
        zsum_weight_patch(jj,jpatch) =  zsum_weight_patch(jj,jpatch)+ zweight(jj,jvegtype)
        IF (pdata(jj,jvegtype).NE.0.) THEN
          zwork(jj,jpatch)= zwork(jj,jpatch) + 1./(log(zdz(jj,jpatch)/ pdata(jj,jvegtype)))**2    &
                            * zweight(jj,jvegtype)
        ENDIF
      ENDDO
    END DO   
!
!-------------------------------------------------------------------------------
!
  CASE default
    CALL abor1_sfx('AV_PGD_PARAM: (1) AVERAGING TYPE NOT ALLOWED')
!
END SELECT
!
!-------------------------------------------------------------------------------
!
!*    4.     End of Averaging
!            ----------------
!
nmask(:,:)=0
jcount(:)=0
DO jpatch=1,ipatch
  DO jj=1,SIZE(zwork,1)
    IF ( zsum_weight_patch(jj,jpatch) >0.) THEN
      jcount(jpatch)=jcount(jpatch)+1
      nmask(jcount(jpatch),jpatch)=jj
    ENDIF
  ENDDO
ENDDO

!*    4.1    Selection of averaging type
!            ---------------------------
!
SELECT CASE (hatype)
!
!-------------------------------------------------------------------------------
!
!*    4.2    Arithmetic averaging
!            --------------------
!
  CASE ('ARI')
!   
    DO jpatch=1,ipatch
!cdir nodep
      DO jj=1,jcount(jpatch)
          ji = nmask(jj,jpatch)
          pfield(ji,jpatch) =  zwork(ji,jpatch) / zsum_weight_patch(ji,jpatch)
      ENDDO
    ENDDO
!
!-------------------------------------------------------------------------------
!
!*    4.3    Inverse averaging
!            -----------------
!
  CASE('INV' )
!
    DO jpatch=1,ipatch
!cdir nodep
      DO jj=1,jcount(jpatch)
        ji = nmask(jj,jpatch)
        pfield(ji,jpatch) = zsum_weight_patch(ji,jpatch) / zwork(ji,jpatch)
      ENDDO
    ENDDO
!-------------------------------------------------------------------------------!
!
!*    4.4    Roughness length averaging
!            --------------------------

!
  CASE('CDN')
!
    DO jpatch=1,ipatch
!cdir nodep
      DO jj=1,jcount(jpatch)
        ji=nmask(jj,jpatch)
        pfield(ji,jpatch) = zdz(ji,jpatch) * exp( - sqrt(zsum_weight_patch(ji,jpatch)/zwork(ji,jpatch)) )
      ENDDO
    ENDDO
!
!-------------------------------------------------------------------------------
!
  CASE default
    CALL abor1_sfx('AV_PGD_PARAM: (2) AVERAGING TYPE NOT ALLOWED')
!
END SELECT
IF (lhook) CALL dr_hook('AV_PGD_PARAM',1,zhook_handle)
!-------------------------------------------------------------------------------
!   
END SUBROUTINE av_pgd_param