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 ( DGU, 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_diag_surf_atm_n, ONLY : diag_surf_atm_t
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_surf_atm_sso_n, ONLY : surf_atm_sso_t
!
USE modd_surf_par,       ONLY : xundef
USE modd_pgd_grid,       ONLY : cgrid, nl, xgrid_par, ngrid_par
USE modd_pgdwork,        ONLY : xsumcover, nsize
USE modd_data_cover_par, ONLY : jpcover, nrock, nsea, nwater, npermsnow
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_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
!
!*    0.1    Declaration of arguments
!            ------------------------
!
!
TYPE(diag_surf_atm_t), INTENT(INOUT) :: dgu
TYPE(data_cover_t), INTENT(INOUT) :: dtco
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
TYPE(surf_atm_t), INTENT(INOUT) :: u
TYPE(surf_atm_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
!
!
INTEGER               :: iluout    ! output listing logical unit
INTEGER               :: iresp     ! Error code after redding
INTEGER               :: jcover    ! loop counter on covers
INTEGER               :: jl        ! 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, DIMENSION(1) :: imaxcover ! index of maximum cover for the given point
INTEGER, DIMENSION(:), POINTER :: imask_cover=>null()
INTEGER, DIMENSION(:), ALLOCATABLE :: imask_sea, imask_water
!
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
!
!---------------------------------------------------------------
!
!*    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 jcover=1,jpcover
      IF (xunif_cover(jcover)/=0.) THEN
        u%LCOVER(jcover) = .true.
        icpt = icpt + 1
        u%XCOVER(:,icpt) = xunif_cover(jcover)
      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, dgu, 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(nl)        )
  ALLOCATE(xsumcover(nl,jpcover))
!
  nsize(:)   = 0.
  xsumcover(:,:) = 0.
  CALL treat_field(ug, u, uss, &
                   hprogram,'SURF  ',yfiletype,'A_COVR',ycover,  &
                     'COVER               '                      ) 

!
!*    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,u%XCOVER(:,:),yfield)
!
!-------------------------------------------------------------------------------
!
!*    5.      Coherence check
!             ---------------
!
  icover = SIZE(u%XCOVER,2)
!
  u%XCOVER(:,:)=u%XCOVER(:,:)/spread(sum(u%XCOVER(:,:),2),2,icover)
!
  DEALLOCATE(nsize    )
  DEALLOCATE(xsumcover)
!
  CALL make_mask_cover(imask_cover,icover)
!
  ALLOCATE(imask_sea(SIZE(nsea)))
  imask_sea(:) = 0
  DO jl=1,SIZE(nsea)
    DO jcover=1,icover
      IF (imask_cover(jcover)==nsea(jl)) imask_sea(jl) = jcover
    ENDDO
  ENDDO
  !
  ALLOCATE(imask_water(SIZE(nwater)))
  imask_water(:) = 0
  DO jl=1,SIZE(nwater)
    DO jcover=1,icover
      IF (imask_cover(jcover)==nwater(jl)) imask_water(jl) = jcover
    ENDDO
  ENDDO
  !
  ipermsnow=0
  DO jcover=1,icover
    IF (imask_cover(jcover)==npermsnow) ipermsnow = jcover
  ENDDO
  !
  ieco2 = 0
  DO jcover=1,icover
    IF (imask_cover(jcover)>300) THEN
      ieco2 = jcover
      EXIT
    ENDIF
  ENDDO
!
!-------------------------------------------------------------------------------
!
!*    6.      Special treatments asked by user
!             --------------------------------
!
! * removes cover with very small coverage
  DO jl=1,SIZE(u%XCOVER,1)
    imaxcover(:) = maxloc(u%XCOVER(jl,:))
    DO jcover=1,icover
      IF (u%XCOVER(jl,jcover)/=0.) THEN
        IF (u%XCOVER(jl,jcover)<=xrm_cover .AND. jcover /= imaxcover(1)) THEN
          u%XCOVER(jl,jcover) = 0.
        END IF
      ENDIF
    END DO
  END DO
  !
  ! * 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,:))
       IF(imask_water(2)/=imaxcover(1).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(u%XCOVER(:,imask_water(jl))<=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(u%XCOVER(:,imask_sea(jl))<=xrm_sea)
           u%XCOVER(:,imask_sea(jl)) = 0.
         ENDWHERE
       ENDIF
     ENDDO          
  ENDIF
  !
  !
  ! * removes cover; replace by sea or inland water if sea > XRM_COAST
  DO jcover=1,icover
    !
    DO jl=1,SIZE(nsea)
      IF (imask_sea(jl)/=0) THEN
        WHERE(u%XCOVER(:,imask_sea(jl))>=xrm_coast)
          u%XCOVER(:,jcover) = 0.
          u%XCOVER(:,imask_sea(jl)) = 1.
        END WHERE 
      ENDIF
    ENDDO
    !
    DO jl=1,SIZE(nwater)
      IF (imask_water(jl)/=0) THEN
        WHERE(u%XCOVER(:,imask_water(jl))>=xrm_coast)
          u%XCOVER(:,jcover) = 0.
          u%XCOVER(:,imask_water(jl)) = 1.
        END WHERE
      ENDIF
    ENDDO 
    !    
  ENDDO
!
!
! * 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
!             ---------------------
!
  ALLOCATE(zcover(nl,icover))
  zcover(:,:) = 0.
!
  icover_old = icover
  icover = 0
!
  ieco2 = 0
!
  u%LCOVER(:) = .false.
  DO jcover=1,icover_old
    icoversum = sum_on_all_procs(hprogram,cgrid,u%XCOVER(:,jcover)/=0., 'COV')
    IF (icoversum>0) THEN
      u%LCOVER(imask_cover(jcover))=.true.
      icover = icover+1
      zcover(:,icover) = u%XCOVER(:,jcover)
      IF (imask_cover(jcover)>300) ieco2 = icover
    ENDIF
  END DO
!
  DEALLOCATE(u%XCOVER)
  ALLOCATE(u%XCOVER(nl,icover))
  u%XCOVER(:,:) = zcover(:,1:icover)
!
  DEALLOCATE(zcover)
  DEALLOCATE(imask_cover)
!
!-------------------------------------------------------------------------------
END IF
!
DEALLOCATE(xunif_cover)
!-------------------------------------------------------------------------------
!
!
IF(.NOT.limp_cover)THEN
        
!*    8.      List of cover present
!             ---------------------
!
  IF (ieco2/=0) THEN
    IF ( sum_on_all_procs(hprogram,cgrid,any(u%XCOVER(:,ieco2:)>0.,dim=2),'COV' ) >0 ) &
      CALL pgd_ecoclimap2_data(dtco, &
                               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))
  !
  ALLOCATE(zdef(icover))
  !
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  WRITE(iluout,fmt=*) &
  '*  Coherence computation between covers and imposed nature fraction *'
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  nsize(:) = 1
  WHERE (u%XNATURE(:).NE.0. .AND. znature(:).EQ.0.) nsize(:)=0
          
  DO jl=1,SIZE(u%XCOVER,1)
    IF (u%XNATURE(jl).EQ.0.) nsize(jl)=-1
  ENDDO
  zdef(:)=0.
  DO jcover=1,icover
    IF (xdata_nature(imask_cover(jcover))/=0.) THEN
      zdef(jcover) = 1.
      EXIT
    ENDIF
  ENDDO
  CALL interpol_field2d(ug, u, &
                        hprogram,iluout,nsize,zcover_nature(:,:),yfield,zdef)  
!
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  WRITE(iluout,fmt=*) &
  '*  Coherence computation between covers and imposed town   fraction *'
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  nsize(:) = 1
  WHERE (u%XTOWN(:).NE.0. .AND. ztown(:).EQ.0.) nsize(:)=0
  DO jl=1,SIZE(u%XCOVER,1)
    IF (u%XTOWN(jl).EQ.0.) nsize(jl)=-1
  ENDDO
  zdef(:)=0.
  DO jcover=1,icover
    IF (xdata_town(imask_cover(jcover))/=0.) THEN
      zdef(jcover) = 1.
      EXIT
    ENDIF
  ENDDO  
  CALL interpol_field2d(ug, u, &
                        hprogram,iluout,nsize,zcover_town(:,:),yfield,zdef) 

  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  WRITE(iluout,fmt=*) &
  '*  Coherence computation between covers and imposed water  fraction *'
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  nsize(:) = 1
  WHERE (u%XWATER(:).NE.0. .AND. zwater(:).EQ.0.) nsize(:)=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
     ELSEIF(u%XWATER(jl)==0.0)THEN
        nsize(jl)=-1
     ENDIF
  ENDDO
  zdef(:)=0.
  DO jcover=1,icover
    IF (xdata_water(imask_cover(jcover))/=0.) THEN
      zdef(jcover) = 1.
      EXIT
    ENDIF
  ENDDO    
  CALL interpol_field2d(ug, u, &
                        hprogram,iluout,nsize,zcover_water(:,:),yfield,pdef=zdef)
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  WRITE(iluout,fmt=*) &
  '*  Coherence computation between covers and imposed sea    fraction *'
  WRITE(iluout,fmt=*) &
  '*********************************************************************'
  nsize(:) = 1
  WHERE (u%XSEA(:).NE.0. .AND. zsea(:).EQ.0.) nsize(:)=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
     ELSEIF(u%XSEA(jl)==0.0)THEN
        nsize(jl)=-1
     ENDIF
  ENDDO
  zdef(:)=0.
  DO jcover=1,icover
    IF (xdata_sea(imask_cover(jcover))/=0.) THEN
      zdef(jcover) = 1.
      EXIT
    ENDIF
  ENDDO    
  CALL interpol_field2d(ug, u, &
                        hprogram,iluout,nsize,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
!
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
!
gpresent = .false.
DO jcover=1,kcover
  IF (pdata_surf(imask_cover(jcover))/=0.) THEN
    gpresent = .true.
    EXIT
  ENDIF
ENDDO
!
IF (any(psurf(:)/=0.)) THEN
  !
  IF (gpresent) THEN
    !
    DO jcover=1,kcover
      IF (imask_cover(jcover)==ksurf) THEN
        kc_surf = jcover
        EXIT
      ENDIF
    ENDDO
    !
  ELSE
    !
    u%LCOVER(ksurf) = .true.
    kcover = kcover + 1
    ALLOCATE(zcover(nl,kcover))
    DO jcover = 1,kcover
      IF (jcover<kcover) THEN
        IF (imask_cover(jcover)<ksurf) cycle
      ENDIF
      kc_surf = jcover
      IF (jcover>1) zcover(:,1:jcover-1) = u%XCOVER(:,1:jcover-1)
      zcover(:,jcover) = 0.
      IF (jcover<kcover) zcover(:,jcover+1:kcover) = u%XCOVER(:,jcover: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
!
END SUBROUTINE fit_covers
!
!------------------------------------------------------
!
SUBROUTINE make_mask_cover(KMASK_COVER,KCOVER)
!
INTEGER, DIMENSION(:), POINTER :: kmask_cover
INTEGER, INTENT(IN) :: kcover
!
INTEGER :: icpt
!
IF (ASSOCIATED(kmask_cover)) DEALLOCATE(kmask_cover)
ALLOCATE(kmask_cover(kcover))
icpt = 0
DO jcover=1,jpcover
  IF (u%LCOVER(jcover)) THEN
    icpt = icpt + 1
    kmask_cover(icpt) = jcover
  ENDIF
ENDDO
!
END SUBROUTINE make_mask_cover
!
END SUBROUTINE pgd_cover