pgd_fieldin.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 pgd_fieldin (DTCO, UG, U, USS, &
                            HPROGRAM,HFIELD,HAREA,HFILE,HFILETYPE,PUNIF,PFIELD,&
                            OPRESENT,PVEGTYPE)
!     ##############################################################
!
!!**** *PGD_FIELDIN* monitor for averaging and interpolations of ISBA physiographic fields
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!   
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson        Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    10/12/97
!!    09/2010 (E. Kourzeneva):   interpolation of the lake depth 
!!                               is not allowed and not necessary
!!
!!    02/2014 (B. Decharme):     interpolation of the lake depth 
!!                               re-allowed but using the nearest point
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
!
USE modd_surfex_mpi, ONLY : nrank, npio, nproc
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_sso_n, ONLY : sso_t
!
USE modd_pgdwork,        ONLY : xall, nsize_all, catype, nsize, xsumval,   &
                                nvalnbr, nvalcount, xvallist, jpvalmax
USE modd_surf_par,       ONLY : xundef
USE modd_pgd_grid,       ONLY : nl
!
USE modd_data_cover_par, ONLY : ntype, lveg_pres, nvegtype
!
USE modi_get_luout
USE modi_treat_field
USE modi_interpol_field
USE modi_pack_same_rank
USE modi_unpack_same_rank
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
USE modi_get_surf_mask_n
!
USE modi_get_type_dim_n
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(sso_t), INTENT(INOUT) :: USS
!
 CHARACTER(LEN=6),  INTENT(IN) :: HPROGRAM  ! Type of program
 CHARACTER(LEN=*),  INTENT(IN) :: HFIELD    ! field name for prints
 CHARACTER(LEN=3),  INTENT(IN) :: HAREA     ! area where field is defined
!                                           ! 'ALL' : everywhere
!                                           ! 'NAT' : on nature
!                                           ! 'TWN' : on town
!                                           ! 'SEA' : on sea
!                                           ! 'WAT' : on inland waters
 CHARACTER(LEN=28), INTENT(IN) :: HFILE     ! data file name
 CHARACTER(LEN=6),  INTENT(INOUT) :: HFILETYPE ! data file type
REAL,              INTENT(IN) :: PUNIF     ! prescribed uniform value for field
REAL, DIMENSION(:,:),INTENT(OUT):: PFIELD    ! physiographic field
LOGICAL, OPTIONAL, INTENT(OUT) :: OPRESENT
REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PVEGTYPE
!
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
INTEGER                        :: ILU    ! expected physical size of full surface array
INTEGER                        :: ILUOUT ! output listing logical unit
INTEGER, DIMENSION(:), POINTER :: IMASK  ! mask for packing from complete field to nature field
INTEGER                        :: IDIM   !
INTEGER :: JI, ICASE
!
REAL, DIMENSION(:), ALLOCATABLE :: ZVEGTYPE
INTEGER :: JJ, JT, JTN
 CHARACTER(LEN=20)   :: YFIELD
 CHARACTER(LEN=6)    :: YMASK
INTEGER             :: INPTS     ! number of points used for interpolation
REAL, DIMENSION(:,:), ALLOCATABLE :: ZFIELD    ! physiographic field on full grid
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
!*    1.      Initializations
!             ---------------
!
IF (lhook) CALL dr_hook('PGD_FIELDIN',0,zhook_handle)
!
!-------------------------------------------------------------------------------
!
!*    2.      Output listing logical unit
!             ---------------------------
!
 CALL get_luout(hprogram,iluout)
!
IF (len_trim(hfile)/=0 .OR. punif/=xundef) THEN
  !
  IF (PRESENT(opresent)) opresent=.true.
  !
  IF (hfiletype=='DIRTYP') THEN
    ALLOCATE(zfield(nl,sum(ntype)))
  ELSE
    ALLOCATE(zfield(nl,1))
  ENDIF
  !-------------------------------------------------------------------------------
  !
  !*    6.      Mask for the field
  !             ------------------
  !
  ymask = ''
  SELECT CASE (harea)
    CASE ('LAN')
          ymask = 'LAND  '
    CASE ('TWN')
          ymask = 'TOWN  '
    CASE ('BLD')
          ymask = 'TOWN '              
    CASE ('NAT')
          ymask = 'NATURE'
    CASE ('SEA')
          ymask = 'SEA   '
    CASE ('WAT')
          ymask = 'WATER '
    CASE DEFAULT
          ymask = 'FULL  '
  END SELECT

  CALL get_type_dim_n(dtco, u, ymask,idim)
  IF (idim/=SIZE(pfield,1)) THEN
     WRITE(iluout,*)'Wrong dimension of MASK: ',idim,SIZE(pfield,1)
     CALL abor1_sfx('PGD_FIELDIN: WRONG DIMENSION OF MASK')
  ENDIF

  ALLOCATE(imask(idim))
  ilu=0
  CALL get_surf_mask_n(dtco, u, ymask,idim,imask,ilu,iluout)
!
ELSE
  !
  IF (PRESENT(opresent)) THEN
    opresent=.false.
    pfield(:,:) = xundef
    IF (lhook) CALL dr_hook('PGD_FIELDIN',1,zhook_handle)
    RETURN
  ENDIF
  !
  WRITE(iluout,*) ' '
  WRITE(iluout,*) '***********************************************************'
  WRITE(iluout,*) '* Error in PGD field preparation of field : ', hfield
  WRITE(iluout,*) '* There is no prescribed value and no input file          *'
  WRITE(iluout,*) '***********************************************************'
  WRITE(iluout,*) ' '
  CALL abor1_sfx('PGD_FIELDIN: NO PRESCRIBED VALUE NOR INPUT FILE FOR '//hfield)
  !
ENDIF
!
!-------------------------------------------------------------------------------
!
!*    3.      Read from file
!             --------------
!
IF (len_trim(hfile)/=0) THEN
!       
!-------------------------------------------------------------------------------
!
!*    4.      Averages the field
!             ------------------
!
  ALLOCATE(nsize_all(u%NDIM_FULL,1))
!
  nsize_all(:,1) = 0
!
  IF (catype=='MAJ') THEN
    ALLOCATE(nvalnbr(u%NDIM_FULL,1))
    ALLOCATE(nvalcount(u%NDIM_FULL,jpvalmax,1))
    ALLOCATE(xvallist(u%NDIM_FULL,jpvalmax,1))
    nvalnbr   = 0
    nvalcount = 0
    xvallist  = xundef
    inpts     = 1
  ELSE
    ALLOCATE(xall(u%NDIM_FULL,1,1))
    xall(:,:,:) = 0.
    inpts       = 3
  ENDIF
!
  IF(hfield=="water depth") THEN
    inpts = 1
  ENDIF
!
  yfield = '                    '
  yfield = hfield(1:min(len(hfield),20))
!
  zfield(:,:) = xundef
!
  CALL treat_field(ug, u, uss, &
                   hprogram,'SURF  ',hfiletype,'A_MESH',hfile,   &
                   yfield,zfield                   )
!
!-------------------------------------------------------------------------------
!
!*    4.      Mask for the interpolations
!             ---------------------------
!
  DO jt=1,SIZE(nsize,2)

    SELECT CASE (harea)

      CASE ('LAN')
        WHERE ((u%XTOWN(:)+u%XNATURE(:))==0. .AND. nsize(:,jt)==0 ) nsize(:,jt) = -1

      CASE ('TWN')
        WHERE (u%XTOWN  (:)==0. .AND. nsize(:,jt)==0 ) nsize(:,jt) = -1

      CASE ('BLD')
        WHERE (u%XTOWN  (:)==0. .AND. nsize(:,jt)==0 ) nsize(:,jt) = -1 

      CASE ('NAT')
        WHERE (u%XNATURE(:)==0. .AND. nsize(:,jt)==0 ) nsize(:,jt) = -1
        !
        ! for fields calculated by vegtype
        IF (PRESENT(pvegtype) .AND.SIZE(nsize,2)>1) THEN
          ! only for the natural part
          IF (u%LECOSG) THEN
            jtn = jt - sum(ntype(1:2))
            nsize(:,1:sum(ntype(1:2))) = -1
          ELSE
            jtn = jt
          ENDIF
          IF ( jtn <= SIZE(pvegtype,2) ) THEN
            ALLOCATE(zvegtype(u%NSIZE_FULL))
            CALL unpack_same_rank(imask,pvegtype(:,jtn),zvegtype)
            WHERE (zvegtype(:)==0. .AND. nsize(:,jt)==0) nsize(:,jt) = -1
            DEALLOCATE(zvegtype)
          ELSE
            nsize(:,jt) = -1
          ENDIF
        ENDIF

      CASE ('SEA')
        WHERE (u%XSEA   (:)==0. .AND. nsize(:,jt)==0 ) nsize(:,jt) = -1

      CASE ('WAT')
        WHERE (u%XWATER (:)==0. .AND. nsize(:,jt)==0 ) nsize(:,jt) = -1

    END SELECT

  ENDDO
!
!-------------------------------------------------------------------------------
!
!*    5.      Interpolation if some points are not initialized (no data for these points)
!             ------------------------------------------------
!
  DO jt=1,SIZE(nsize,2)

    IF (.NOT.u%LECOSG.AND.jt>nvegtype) EXIT
    
    !multitype input file
    IF (SIZE(zfield,2)>1) THEN

      !fields defined only on the not-bare soil part
      IF ( (yfield(1:3)=='LAI'.OR.yfield(1:10)=='ALBNIR_VEG'.OR.yfield(1:10)=='ALBVIS_VEG') .OR. &
           (SIZE(zfield,2)>1.AND.yfield(1:6)=='H_TREE') ) THEN

        IF ( (.NOT.u%LECOSG.AND.jt<=3).OR.(u%LECOSG.AND.jt<=sum(ntype(1:2))+3) ) THEN
          zfield(:,jt) = 0.
          nsize(:,jt) = 1
        ENDIF

        ! height of trees only defined for tree vegtypes
        IF (yfield(1:6)=='H_TREE') THEN
          IF ((.NOT.u%LECOSG.AND.((jt>=7 .AND. jt<=12) .OR. (jt>=18 .AND. jt<=19))).OR. &
              (     u%LECOSG.AND.((jt<=(sum(ntype(1:2))+3).OR.jt>=(sum(ntype(1:2))+13)).AND.jt/=22) ) ) THEN
            zfield(:,jt) = 0.
            nsize(:,jt) = 1.
          ENDIF
        ENDIF

      ENDIF

      ! if the cover / vegtype is not present on the area
      !ICASE = 0
      !IF (.NOT.U%LECOSG) THEN
      !  IF (.NOT.LVEG_PRES(JT)) ICASE=1
      !ENDIF
      !IF( (U%LECOSG.AND..NOT.U%LCOVER(JT)) .OR. ICASE==1 ) THEN
      !  ZFIELD(:,JT) = 0.
      !  NSIZE (:,JT) = 1.
      !ENDIF

    ENDIF

    IF (punif/=xundef) THEN
      CALL interpol_field(ug, u, hprogram,iluout,nsize(:,jt),zfield(:,jt),hfield,pdef=punif,knpts=inpts)
    ELSE
      CALL interpol_field(ug, u, hprogram,iluout,nsize(:,jt),zfield(:,jt),hfield)
    ENDIF

  ENDDO        
!
  DEALLOCATE(nsize    )
  DEALLOCATE(xsumval  )
!
!-------------------------------------------------------------------------------
!
ELSEIF (punif/=xundef) THEN
!
!*    3.1     Use of the presribed field
!             --------------------------
!
  zfield(:,:) = punif
!
!
END IF
!
! only the case nature is treated for now, to adapt for town later
IF (harea=='NAT'.AND.SIZE(zfield,2)>SIZE(pfield,2).AND.u%LECOSG) THEN
  CALL pack_same_rank(imask,zfield(:,sum(ntype(1:2))+1:sum(ntype(1:3))),pfield(:,:))
ELSE
  CALL pack_same_rank(imask,zfield(:,1:SIZE(pfield,2)),pfield(:,:))
ENDIF
!
DEALLOCATE(zfield)
DEALLOCATE(imask)
!
IF (lhook) CALL dr_hook('PGD_FIELDIN',1,zhook_handle)

!
!-------------------------------------------------------------------------------
!
END SUBROUTINE pgd_fieldin