pgd_cover.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_cover ( DTCO, UG, U, USS, HPROGRAM, ORM_RIVER)
!     ##############################################################
!
!!**** *PGD_COVER* monitor for averaging and interpolations of cover fractions
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!   
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson        Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    10/12/97
!!    B. Decharme  06/2008  limit of coast coverage under which the coast is replaced by sea or inland water
!!    B. Decharme  06/2009  remove lack and sea as the user want
!!    B. Decharme  07/2009  compatibility between Surfex and Orca (Nemo) grid (Earth Model)
!!    B. Decharme  07/2012  if sea or water imposed to 1 in a grid cell: no extrapolation
!!    B. Decharme  02/2014  Add LRM_RIVER and remove lake over antarctica
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
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_surfex_mpi, ONLY : nrank, npio, nproc, ncomm
USE modd_surf_par,       ONLY : xundef
USE modd_pgd_grid,       ONLY : cgrid, nl, xgrid_par
USE modd_pgdwork,        ONLY : xall, nsize_all, nsize, xsumval, xprec
USE modd_data_cover_par, ONLY : jpcover, nrock, nsea, nwater, npermsnow, lveg_pres
USE modd_data_cover,     ONLY : xdata_town, xdata_sea, xdata_nature, xdata_water
!
USE modi_get_luout
USE mode_gridtype_gauss
USE mode_gridtype_lonlat_reg
!
USE mode_read_surf_cov, ONLY : read_surf_cov
!
USE modi_treat_field
USE modi_interpol_field2d
USE modi_convert_cover_frac
!
USE modi_read_lcover
USE modi_sum_on_all_procs
!
USE modi_make_lcover
USE modi_read_nam_pgd_cover
!
USE modi_init_io_surf_n
USE modi_end_io_surf_n
!
USE modi_abor1_sfx
!
USE modi_pgd_ecoclimap2_data
!
#ifdef SFX_ASC
USE modd_io_surf_asc, ONLY : cfilein
#endif
#ifdef SFX_FA
USE modd_io_surf_fa,  ONLY : cfilein_fa
#endif
#ifdef SFX_LFI
USE modd_io_surf_lfi, ONLY : cfilein_lfi
#endif
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include "mpif.h"
#endif
!
!*    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
LOGICAL,             INTENT(OUT)   :: ORM_RIVER    ! delete river coverage (default = false)
!
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
 CHARACTER(LEN=10)       :: YFIELD
 CHARACTER(LEN=28)       :: YCOVER      ! file name for cover types
 CHARACTER(LEN=6)        :: YFILETYPE   ! data file type
!
REAL                     :: XRM_COVER   ! limit of coverage under which the
                                        ! cover is removed. Default is 1.E-6
REAL                     :: XRM_COAST   ! limit of coast coverage under which
                                        ! the coast is replaced by sea. Default is 1.
REAL                     :: XRM_LAKE    ! limit of inland lake coverage under which
                                        ! the water is removed. Default is 0.0                     
REAL                     :: XRM_SEA     ! limit of sea coverage under which
                                        ! the sea is removed. Default is 0.0
REAL                     :: XLAT_ANT    ! Lattitude limit from Orca grid (Antartic)
!
REAL, DIMENSION(:), ALLOCATABLE :: ZDEF
REAL, DIMENSION(:), ALLOCATABLE :: ZLAT
REAL, DIMENSION(:), ALLOCATABLE :: XUNIF_COVER ! value of each cover (cover will be
!                                                uniform on the horizontal)
REAL, DIMENSION(:), ALLOCATABLE :: ZSEA   !to check compatibility between 
REAL, DIMENSION(:), ALLOCATABLE :: ZWATER !prescribed fractions and ECOCLIMAP
REAL, DIMENSION(:), ALLOCATABLE :: ZNATURE
REAL, DIMENSION(:), ALLOCATABLE :: ZTOWN
REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOVER_NATURE, ZCOVER_TOWN, ZCOVER_SEA, ZCOVER_WATER, &
                                     ZCOVER, ZCOVER2
!
LOGICAL, DIMENSION(:), ALLOCATABLE :: GCOVER, GCOVER2
!
INTEGER :: INFOMPI, JPROC
INTEGER               :: ILUOUT    ! output listing logical unit
INTEGER               :: IRESP     ! Error code after redding
INTEGER               :: JCOV    ! loop counter on covers
INTEGER               :: JL, JI    ! loop counter on horizontal points
INTEGER               :: ICOVER, ICOVERSUM, ICOVER_OLD, ICPT  ! 0 if cover is not present, >1 if present somewhere
INTEGER               :: IPERMSNOW, IECO2 
INTEGER               :: IC_NAT, IC_TWN, IC_WAT, IC_SEA
INTEGER :: ICPT1, ICPT2, ICPT_TOT
!
INTEGER :: IMAXCOVER ! index of maximum cover for the given point
INTEGER, DIMENSION(:), POINTER :: IMASK_COVER=>null()
INTEGER, DIMENSION(:), ALLOCATABLE :: IMASK_SEA, IMASK_WATER
!
LOGICAL, DIMENSION(:,:), ALLOCATABLE :: GCOVER_ALL
LOGICAL, DIMENSION(:), ALLOCATABLE :: GCORINE
LOGICAL                  :: LORCA_GRID  ! flag to compatibility between Surfex and Orca grid 
                                        ! (Earth Model over Antarctic)
LOGICAL                  :: LIMP_COVER  ! Imposed values for Cover from another PGD file
!
LOGICAL                  :: GPRESENT
!
LOGICAL                  :: LRM_RIVER   ! delete inland river coverage. Default is false
!
REAL, PARAMETER          :: ZLAT_ANT_WATER = -60. ! Lattitude limit to delete lake over antarctica
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!
!---------------------------------------------------------------
!
!*    1.      Initializations
!             ---------------
!
IF (lhook) CALL dr_hook('PGD_COVER',0,zhook_handle)
!
 CALL get_luout(hprogram,iluout)
!
ALLOCATE(u%LCOVER   (jpcover))
ALLOCATE(xunif_cover(jpcover))
!
u%LCOVER    = .false.
xunif_cover = xundef
!
ieco2 = 0
!
!-------------------------------------------------------------------------------
!
!*    2.      Input file for cover types
!             --------------------------
!
 CALL read_nam_pgd_cover(hprogram, ycover, yfiletype, xunif_cover,  &
                         xrm_cover, xrm_coast, xrm_lake, lrm_river, &
                         xrm_sea, lorca_grid, xlat_ant, limp_cover )  
!
!-------------------------------------------------------------------------------
!
!*    3.      Uniform field is prescribed
!             ---------------------------
!-------------------------------------------------------------------------------
!
IF (any(xunif_cover/=0.)) THEN
!
!*    3.1     Verification of the total input cover fractions
!             -----------------------------------------------
!
  IF (abs(sum(xunif_cover)-1.)>1.e-6) THEN
    WRITE(iluout,*) ' '
    WRITE(iluout,*) '***************************************************'
    WRITE(iluout,*) '* Error in COVER fractions preparation            *'
    WRITE(iluout,*) '* The prescribed covers does not fit              *'
    WRITE(iluout,*) '* The sum of all cover must be equal to 1 exactly *'
    WRITE(iluout,*) '***************************************************'
    WRITE(iluout,*) ' '
    CALL abor1_sfx('PGD_COVER: SUM OF ALL COVER FRACTIONS MUST BE 1.')
!
!*    3.2     Use of the presribed cover fractions
!             ------------------------------------
!
  ELSE
    icover = count(xunif_cover(:)/=0.)
    ALLOCATE(u%XCOVER(nl,icover))
    icpt = 0
    DO jcov=1,jpcover
      IF (xunif_cover(jcov)/=0.) THEN
        u%LCOVER(jcov) = .true.
        icpt = icpt + 1
        u%XCOVER(:,icpt) = xunif_cover(jcov)
      ENDIF
    END DO
    u%XCOVER(:,:) = u%XCOVER(:,:) / spread(sum(u%XCOVER(:,:),2),2,icover)
  END IF
!
!*    3.3     No data
!             -------
!
ELSEIF (len_trim(ycover)==0) THEN
  WRITE(iluout,*) ' '
  WRITE(iluout,*) '***********************************************************'
  WRITE(iluout,*) '* Error in COVER fractions preparation                    *'
  WRITE(iluout,*) '* There is no prescribed cover fraction and no input file *'
  WRITE(iluout,*) '***********************************************************'
  WRITE(iluout,*) ' '
  CALL abor1_sfx('PGD_COVER: NO PRESCRIBED COVER NOR INPUT FILE')
!
!-------------------------------------------------------------------------------
ELSEIF(limp_cover)THEN !LIMP_COVER (impose cover from input file at the same resolution)
!
  IF(yfiletype=='NETCDF')THEN
    CALL abor1_sfx('Use another format than netcdf for cover input file with LIMP_COVER')
  ELSE
#ifdef SFX_ASC
    cfilein     = adjustl(adjustr(ycover)//'.txt')
#endif
#ifdef SFX_FA
    cfilein_fa  = adjustl(adjustr(ycover)//'.fa')
#endif
#ifdef SFX_LFI
    cfilein_lfi = adjustl(ycover)
#endif
 CALL init_io_surf_n(dtco, u, &
                        yfiletype,'FULL  ','SURF  ','READ ')
  ENDIF
!
  ALLOCATE(u%LCOVER(jpcover))
  CALL read_lcover(yfiletype,u%LCOVER)
!
  CALL read_surf_cov(yfiletype,'COVER',u%XCOVER(:,:),u%LCOVER,iresp)
!
  CALL end_io_surf_n(yfiletype)
!
ELSE 
!-------------------------------------------------------------------------------
!
!*    3.      Averages the field
!             ------------------
!
  ALLOCATE(nsize_all(u%NDIM_FULL,1)      )
  ALLOCATE(xall(u%NDIM_FULL,1,2)    )
!
  nsize_all(:,:) = 0
  xall(:,:,:) = 0.
  CALL treat_field(ug, u, uss, &
                   hprogram,'SURF  ',yfiletype,'A_COVR',ycover, 'COVER               ' ) 
!
  DEALLOCATE(xsumval  )
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!*    4.      Interpolation if some points are not initialized (no data for these points) (same time)
!             ---------------------------------------------------------------------------------------
!
  WRITE(yfield,fmt='(A)') 'covers'
  CALL interpol_field2d(ug, u, hprogram,iluout,nsize(:,1), u%XCOVER(:,:),yfield)
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!*    5.      Coherence check
!             ---------------
!
  icover = SIZE(u%XCOVER,2)
!
  u%XCOVER(:,:) = u%XCOVER(:,:) / spread(sum(u%XCOVER(:,:),2),2,icover)
!
  DEALLOCATE(nsize    )
!
  CALL make_mask_cover(imask_cover,icover)
!
  ALLOCATE(imask_sea(SIZE(nsea)))
  imask_sea(:) = 0
  DO jl=1,SIZE(nsea)
    DO jcov=1,icover
      IF (imask_cover(jcov)==nsea(jl)) imask_sea(jl) = jcov
    ENDDO
  ENDDO
  !
  ALLOCATE(imask_water(SIZE(nwater)))
  imask_water(:) = 0
  DO jl=1,SIZE(nwater)
    DO jcov=1,icover
      IF (imask_cover(jcov)==nwater(jl)) imask_water(jl) = jcov
    ENDDO
  ENDDO
  !
  ipermsnow=0
  DO jcov=1,icover
    IF (imask_cover(jcov)==npermsnow) ipermsnow = jcov
  ENDDO
  !
  ieco2 = 0
  DO jcov=1,icover
    IF (imask_cover(jcov)>300) THEN
      ieco2 = jcov
      EXIT
    ENDIF
  ENDDO
!
!
!-------------------------------------------------------------------------------
!
!*    6.      Special treatments asked by user
!             --------------------------------
!
 CALL get_rmcov_omp(icover,u%XCOVER)
  !
  ! * removes River if the user want
  orm_river=lrm_river
  IF(lrm_river.AND.imask_water(2)/=0)THEN
    DO jl=1,SIZE(u%XCOVER,1)
       imaxcover = maxloc(u%XCOVER(jl,:),1)
       IF(imask_water(2)/=imaxcover.AND.u%XCOVER(jl,imask_water(2))>0.)THEN
         u%XCOVER(jl,imask_water(2)) = 0.
       ENDIF
    ENDDO
  ENDIF
  !
  ! * removes lake as the user want
  IF(xrm_lake>0.0)THEN
     DO jl=1,SIZE(nwater)
       IF (imask_water(jl)/=0) THEN
         WHERE(anint(u%XCOVER(:,imask_water(jl))*xprec)/xprec<=xrm_lake)
           u%XCOVER(:,imask_water(jl)) = 0.
         ENDWHERE
       ENDIF
     ENDDO          
  ENDIF
  !
  ! * removes sea as the user want
  IF(xrm_sea>0.0)THEN
     DO jl=1,SIZE(nsea)
       IF (imask_sea(jl)/=0) THEN
         WHERE(anint(u%XCOVER(:,imask_sea(jl))*xprec)/xprec<=xrm_sea)
           u%XCOVER(:,imask_sea(jl)) = 0.
         ENDWHERE
       ENDIF
     ENDDO          
  ENDIF
  !
  !
  ! * removes cover; replace by sea or inland water if sea > XRM_COAST
  IF (xrm_coast<1.) THEN
    !
    DO jl=1,SIZE(nsea)
      IF (imask_sea(jl)/=0) THEN
        DO ji=1,SIZE(u%XCOVER,1)   
          IF (anint(u%XCOVER(ji,imask_sea(jl))*xprec)/xprec>=xrm_coast .AND. &
                u%XCOVER(ji,imask_sea(jl))/=1.) THEN
            u%XCOVER(ji,:) = 0.
            u%XCOVER(ji,imask_sea(jl)) = 1.
          ENDIF
        ENDDO
      ENDIF
    ENDDO
    !    
    DO jl=1,SIZE(nwater)
      IF (imask_water(jl)/=0) THEN
         DO ji=1,SIZE(u%XCOVER,1)
          IF (anint(u%XCOVER(ji,imask_water(jl))*xprec)/xprec>=xrm_coast .AND. &
                u%XCOVER(ji,imask_water(jl))/=1. ) THEN
            u%XCOVER(ji,:) = 0.
            u%XCOVER(ji,imask_water(jl)) = 1.
          ENDIF
        ENDDO
      ENDIF
    ENDDO     
    !    
  ENDIF
!
!
! * Compatibility between Surfex and Orca grid 
!   (Earth Model over water bodies and Antarctic)
!
  IF(lorca_grid.AND.(cgrid=='GAUSS     '.OR.cgrid=='LONLAT REG'))THEN
!
    ALLOCATE(zlat(nl))
    IF (cgrid=='GAUSS     ') CALL get_gridtype_gauss(xgrid_par,plat=zlat)
    IF (cgrid=='LONLAT REG') CALL get_gridtype_lonlat_reg(xgrid_par,plat=zlat)
!
      DO jl=1,SIZE(nsea)
        IF (imask_sea(jl)/=0.AND.ipermsnow/=0) THEN
          WHERE(zlat(:)<xlat_ant.AND.u%XCOVER(:,imask_sea(jl))>0.0)
            u%XCOVER(:,ipermsnow) = 1.0
            u%XCOVER(:,imask_sea(jl))  = 0.0
          ENDWHERE 
        ENDIF
      ENDDO

      DO jl=1,SIZE(nwater)
        IF (imask_water(jl)/=0.AND.ipermsnow/=0) THEN
          WHERE(zlat(:)<zlat_ant_water.AND.u%XCOVER(:,imask_water(jl))>0.0)
            u%XCOVER(:,ipermsnow)  = 1.0
            u%XCOVER(:,imask_water(jl)) = 0.0
          ENDWHERE
        ENDIF
      ENDDO
!  
    DEALLOCATE(zlat)
!
  ENDIF
!
!-------------------------------------------------------------------------------
!
!*    7.      Coherence check
!             ---------------
!
  u%XCOVER(:,:)=u%XCOVER(:,:)/spread(sum(u%XCOVER(:,:),2),2,icover)
!
  DEALLOCATE(imask_sea)
  DEALLOCATE(imask_water)
!
!*    8.      List of cover present
!             ---------------------
!
  u%LCOVER(:) = .false.
  DO jcov=1,icover
    IF (any(u%XCOVER(:,jcov)/=0.)) u%LCOVER(imask_cover(jcov)) = .true.
  ENDDO
!
 CALL make_lcover(u%LCOVER)
!
  icover_old = icover
  icover = count(u%LCOVER)
!
  IF (hprogram=='MESONH'.OR.icover<icover_old) THEN

    IF (icover/=icover_old) THEN
      ALLOCATE(zcover(nl,icover_old))
      zcover(:,:) = u%XCOVER(:,:)
      DEALLOCATE(u%XCOVER)
      ALLOCATE(u%XCOVER(nl,icover))
      u%XCOVER(:,:) = 0.
      icpt = 0
      DO jcov=1,icover_old
        IF (u%LCOVER(imask_cover(jcov))) THEN
          icpt = icpt + 1
          u%XCOVER(:,icpt) = zcover(:,jcov)
        ENDIF
      ENDDO
      DEALLOCATE(zcover)
    ENDIF

    IF (hprogram=='MESONH') THEN

      DO jcov=1,icover
        icoversum = sum_on_all_procs(hprogram,cgrid,u%XCOVER(:,jcov)/=0., 'COV')
        IF (icoversum==0) u%LCOVER(imask_cover(jcov) )= .false.
      ENDDO

      icover_old = icover
      icover = count(u%LCOVER)
      IF (icover/=icover_old) THEN
        ALLOCATE(zcover(nl,icover_old))
        zcover(:,:) = u%XCOVER(:,:)
        DEALLOCATE(u%XCOVER)
        ALLOCATE(u%XCOVER(nl,icover))
        u%XCOVER(:,:) = 0.
        icpt = 0
        DO jcov=1,icover_old
          IF (u%LCOVER(imask_cover(jcov))) THEN
            icpt = icpt + 1
            u%XCOVER(:,icpt) = zcover(:,jcov)
          ENDIF
        ENDDO
        DEALLOCATE(zcover)
      ENDIF

    ENDIF

  ENDIF
!
  ieco2 = 0
  IF (any(u%LCOVER(301:))) ieco2=1
!
  DEALLOCATE(imask_cover)
!
!-------------------------------------------------------------------------------
END IF
!
DEALLOCATE(xunif_cover)
!-------------------------------------------------------------------------------
!
!
IF(.NOT.limp_cover)THEN
        
!*    8.      List of cover present
!             ---------------------
!
  IF (ieco2/=0) THEN
      CALL pgd_ecoclimap2_data(dtco%NYEAR, dtco%XDATA_VEGTYPE, hprogram)
  ENDIF
!
!-------------------------------------------------------------------------------
ENDIF
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!*    9.      Land - sea fractions
!             --------------------
!
IF (.NOT.ASSOCIATED(u%XSEA)) THEN

  ALLOCATE(u%XSEA   (nl))
  ALLOCATE(u%XWATER (nl))
  ALLOCATE(u%XNATURE(nl))
  ALLOCATE(u%XTOWN  (nl))
  CALL convert_cover_frac(dtco, u%XCOVER,u%LCOVER,u%XSEA,u%XNATURE,u%XTOWN,u%XWATER)

ELSE
  !
  icover = SIZE(u%XCOVER,2)
  !
  CALL make_mask_cover(imask_cover,icover)
  !
!if fractions are prescribed, it has to be verified that the locations of
!ECOCLIMAP covers are compatible with the fractions of surface types
  ALLOCATE(zsea(nl))
  ALLOCATE(zwater(nl))
  ALLOCATE(znature(nl))
  ALLOCATE(ztown(nl))
  CALL convert_cover_frac(dtco, u%XCOVER,u%LCOVER,zsea,znature,ztown,zwater)
  !
  CALL fit_covers(xdata_nature,u%XNATURE,4,icover,ic_nat)
  CALL fit_covers(xdata_town  ,u%XTOWN  ,7,icover,ic_twn)
  CALL fit_covers(xdata_water ,u%XWATER ,2,icover,ic_wat)
  CALL fit_covers(xdata_sea   ,u%XSEA   ,1,icover,ic_sea)
  !
  ALLOCATE(zcover_nature(nl,icover))
  ALLOCATE(zcover_town(nl,icover))
  ALLOCATE(zcover_sea(nl,icover))
  ALLOCATE(zcover_water(nl,icover))
  !
  zcover_nature(:,:) = u%XCOVER(:,:)
  zcover_town(:,:) = u%XCOVER(:,:)
  zcover_sea(:,:) = u%XCOVER(:,:)
  zcover_water(:,:) = u%XCOVER(:,:)
  !
  ALLOCATE(nsize(nl,1))
  !
  ALLOCATE(zdef(icover))
  !
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  WRITE(iluout,fmt=*) &
  '*  Coherence computation between covers and imposed nature fraction *'
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  nsize(:,1) = 1
  WHERE (u%XNATURE(:).NE.0. .AND. znature(:).EQ.0.) nsize(:,1)=0
          
  DO jl=1,SIZE(u%XCOVER,1)
    IF (u%XNATURE(jl).EQ.0.) nsize(jl,1)=-1
  ENDDO
  zdef(:)=0.
  DO jcov=1,icover
    IF (xdata_nature(imask_cover(jcov))/=0.) THEN
      zdef(jcov) = 1.
      EXIT
    ENDIF
  ENDDO
  CALL interpol_field2d(ug, u, &
                        hprogram,iluout,nsize(:,1),zcover_nature(:,:),yfield,zdef)  
!
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  WRITE(iluout,fmt=*) &
  '*  Coherence computation between covers and imposed town   fraction *'
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  nsize(:,1) = 1
  WHERE (u%XTOWN(:).NE.0. .AND. ztown(:).EQ.0.) nsize(:,1)=0
  DO jl=1,SIZE(u%XCOVER,1)
    IF (u%XTOWN(jl).EQ.0.) nsize(jl,1)=-1
  ENDDO
  zdef(:)=0.
  DO jcov=1,icover
    IF (xdata_town(imask_cover(jcov))/=0.) THEN
      zdef(jcov) = 1.
      EXIT
    ENDIF
  ENDDO  
  CALL interpol_field2d(ug, u, &
                        hprogram,iluout,nsize(:,1),zcover_town(:,:),yfield,zdef) 

  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  WRITE(iluout,fmt=*) &
  '*  Coherence computation between covers and imposed water  fraction *'
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  nsize(:,1) = 1
  WHERE (u%XWATER(:).NE.0. .AND. zwater(:).EQ.0.) nsize(:,1)=0
! if water imposed to 1 in a grid cell: no extrapolation
  DO jl=1,SIZE(u%XCOVER,1)
     IF(u%XWATER(jl)==1.0)THEN
        zcover_water(jl,:)=0.0             
        zcover_water(jl,ic_wat)=1.0
        nsize(jl,1)=1
     ELSEIF(u%XWATER(jl)==0.0)THEN
        nsize(jl,1)=-1
     ENDIF
  ENDDO
  zdef(:)=0.
  DO jcov=1,icover
    IF (xdata_water(imask_cover(jcov))/=0.) THEN
      zdef(jcov) = 1.
      EXIT
    ENDIF
  ENDDO    
  CALL interpol_field2d(ug, u, &
                        hprogram,iluout,nsize(:,1),zcover_water(:,:),yfield,pdef=zdef)
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  WRITE(iluout,fmt=*) &
  '*  Coherence computation between covers and imposed sea    fraction *'
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  nsize(:,1) = 1
  WHERE (u%XSEA(:).NE.0. .AND. zsea(:).EQ.0.) nsize(:,1)=0
! if sea imposed to 1 in a grid cell: no extrapolation          
  DO jl=1,SIZE(u%XCOVER,1)
     IF(u%XSEA(jl)==1.0)THEN
        zcover_sea(jl,:)=0.0             
        zcover_sea(jl,ic_sea)=1.0
        nsize(jl,1)=1
     ELSEIF(u%XSEA(jl)==0.0)THEN
        nsize(jl,1)=-1
     ENDIF
  ENDDO
  zdef(:)=0.
  DO jcov=1,icover
    IF (xdata_sea(imask_cover(jcov))/=0.) THEN
      zdef(jcov) = 1.
      EXIT
    ENDIF
  ENDDO    
  CALL interpol_field2d(ug, u, &
                        hprogram,iluout,nsize(:,1),zcover_sea(:,:),yfield,pdef=zdef)
  !
  u%XCOVER(:,:) = u%XCOVER(:,:) + 0.001 * ( zcover_nature(:,:) + zcover_town(:,:) + &
                                            zcover_water(:,:) + zcover_sea(:,:) )
  !
  u%XCOVER(:,:)=u%XCOVER(:,:)/spread(sum(u%XCOVER(:,:),2),2,icover)
  !
  DEALLOCATE(zcover_nature)
  DEALLOCATE(zcover_town  )
  DEALLOCATE(zcover_water )
  DEALLOCATE(zcover_sea   )
  !
  DEALLOCATE(nsize    )
  DEALLOCATE(zsea     )
  DEALLOCATE(zwater   )
  DEALLOCATE(znature  )
  DEALLOCATE(ztown    )
  !
  DEALLOCATE(zdef)
  DEALLOCATE(imask_cover)
  !
ENDIF
!
lveg_pres(:) = .false.
IF (.NOT.u%LECOSG) THEN
  DO jcov = 1,jpcover
    IF (u%LCOVER(jcov)) THEN
      WHERE(dtco%XDATA_VEGTYPE(jcov,:) > 0.) lveg_pres(1:SIZE(dtco%XDATA_VEGTYPE,2)) = .true.
    ENDIF
  ENDDO
ENDIF
!
u%NSIZE_NATURE    = count(u%XNATURE(:) > 0.0)
u%NSIZE_WATER     = count(u%XWATER (:) > 0.0)
u%NSIZE_SEA       = count(u%XSEA   (:) > 0.0)
u%NSIZE_TOWN      = count(u%XTOWN  (:) > 0.0)
u%NSIZE_FULL      = nl
!
u%NDIM_NATURE    = sum_on_all_procs(hprogram,cgrid,u%XNATURE(:) > 0., 'DIM')
u%NDIM_WATER     = sum_on_all_procs(hprogram,cgrid,u%XWATER (:) > 0., 'DIM')
u%NDIM_SEA       = sum_on_all_procs(hprogram,cgrid,u%XSEA   (:) > 0., 'DIM')
u%NDIM_TOWN      = sum_on_all_procs(hprogram,cgrid,u%XTOWN  (:) > 0., 'DIM')
!
IF (lhook) CALL dr_hook('PGD_COVER',1,zhook_handle)
!-------------------------------------------------------------------------------
CONTAINS
!
SUBROUTINE fit_covers(PDATA_SURF,PSURF,KSURF,KCOVER,KC_SURF)
!
REAL, DIMENSION(:), INTENT(IN) :: PDATA_SURF
REAL, DIMENSION(:), INTENT(IN) :: PSURF
INTEGER, INTENT(IN) :: KSURF
INTEGER, INTENT(INOUT) :: KCOVER
INTEGER, INTENT(OUT) :: KC_SURF
!
LOGICAL :: GPRESENT
REAL :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('PGD_COVER:FIT_COVERS',0,zhook_handle)
!
gpresent = .false.
DO jcov=1,kcover
  IF (pdata_surf(imask_cover(jcov))/=0.) THEN
    gpresent = .true.
    EXIT
  ENDIF
ENDDO
!
IF (any(psurf(:)/=0.)) THEN
  !
  IF (gpresent) THEN
    !
    DO jcov=1,kcover
      IF (imask_cover(jcov)==ksurf) THEN
        kc_surf = jcov
        EXIT
      ENDIF
    ENDDO
    !
  ELSE
    !
    u%LCOVER(ksurf) = .true.
    kcover = kcover + 1
    ALLOCATE(zcover(nl,kcover))
    DO jcov = 1,kcover
      IF (jcov<kcover) THEN
        IF (imask_cover(jcov)<ksurf) cycle
      ENDIF
      kc_surf = jcov
      IF (jcov>1) zcover(:,1:jcov-1) = u%XCOVER(:,1:jcov-1)
      zcover(:,jcov) = 0.
      IF (jcov<kcover) zcover(:,jcov+1:kcover) = u%XCOVER(:,jcov:kcover-1)
      EXIT
    ENDDO
    DEALLOCATE(u%XCOVER)
    ALLOCATE(u%XCOVER(nl,kcover))
    u%XCOVER(:,:) = zcover(:,:)
    DEALLOCATE(zcover)
    !
    CALL make_mask_cover(imask_cover,kcover)
    !
  ENDIF
  !
ENDIF
!
IF (lhook) CALL dr_hook('PGD_COVER:FIT_COVERS',1,zhook_handle)
!
END SUBROUTINE fit_covers
!
!------------------------------------------------------
!
SUBROUTINE make_mask_cover(KMASK_COVER,KCOVER)
!
INTEGER, DIMENSION(:), POINTER :: KMASK_COVER
INTEGER, INTENT(IN) :: KCOVER
!
INTEGER :: ICPT
REAL :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('PGD_COVER:MAKE_MASK_COVER',0,zhook_handle)
!
IF (ASSOCIATED(kmask_cover)) DEALLOCATE(kmask_cover)
ALLOCATE(kmask_cover(kcover))
icpt = 0
DO jcov=1,jpcover
  IF (u%LCOVER(jcov)) THEN
    icpt = icpt + 1
    kmask_cover(icpt) = jcov
  ENDIF
ENDDO
!
IF (lhook) CALL dr_hook('PGD_COVER:MAKE_MASK_COVER',1,zhook_handle)
!
END SUBROUTINE make_mask_cover
!
!------------------------------------------------------
!
SUBROUTINE get_rmcov_omp(KCOVER,PCOVER)
!
USE modd_surfex_omp, ONLY : nblocktot
!
USE modi_gather_and_write_mpi
USE modi_read_and_send_mpi
!
INTEGER, INTENT(IN) :: KCOVER
REAL, DIMENSION(:,:), INTENT(INOUT) :: PCOVER
!
REAL, DIMENSION(U%NDIM_FULL,SIZE(PCOVER,2)) :: ZCOVER_ALL
INTEGER, DIMENSION(U%NDIM_FULL) :: IMAXCOVER_ALL
INTEGER, DIMENSION(U%NSIZE_FULL) :: IMAXCOVER
INTEGER :: JK, JCOV, ISIZE_OMP
REAL :: ZHOOK_HANDLE_OMP
!
isize_omp = max(1,SIZE(pcover,1)/nblocktot)
!
 CALL gather_and_write_mpi(pcover,zcover_all)
IF (nrank==npio) THEN
  DO jl = 1,u%NDIM_FULL
    imaxcover_all(jl) = maxloc(zcover_all(jl,:),1)
  ENDDO
ENDIF
!
 CALL read_and_send_mpi(imaxcover_all,imaxcover)
!
! * removes cover with very small coverage
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (lhook) CALL dr_hook('PGD_COVER:GET_RMCOV_OMP',0,zhook_handle_omp)
!$OMP DO SCHEDULE(STATIC,ISIZE_OMP) PRIVATE(JL,JCOV)
DO jl=1,SIZE(pcover,1)
  DO jcov=1,kcover
    IF (jcov /= imaxcover(jl)) THEN
      IF (anint(pcover(jl,jcov)*xprec)/xprec<=xrm_cover ) pcover(jl,jcov) = 0.
    ENDIF
  END DO
END DO
!$OMP END DO
IF (lhook) CALL dr_hook('PGD_COVER:GET_RMCOV_OMP',1,zhook_handle_omp)
!$OMP END PARALLEL
!
END SUBROUTINE get_rmcov_omp
!
END SUBROUTINE pgd_cover