prep_hor_isba_field.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 prep_hor_isba_field (DTCO, IG, I, UG, U, USS, &
                                hprogram,hsurf,hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,okey)
!     #################################################################################
!
!!****  *PREP_HOR_ISBA_FIELD* - reads, interpolates and prepares an ISBA field
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      P. Le Moigne 10/2005, Phasage Arome
!!      P. Le Moigne 03/2007, Ajout initialisation par ascllv
!!      B. Decharme  01/2009, Optional Arpege deep soil temperature initialization
!!      M. Lafaysse  07/2012, allow netcdf input files
!!      B. Decharme  07/2012, Bug init uniform snow
!!      M. Lafaysse 11/2012,  snow liquid water content
!!      B. Decharme  03/2014, external init with FA files
!!                            new vertical interpolation
!!      P Samuelsson 10/2014, MEB
!!------------------------------------------------------------------
!
!
!
!
!
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_isba_grid_n, ONLY : isba_grid_t
USE modd_isba_n, ONLY : isba_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_prep,     ONLY : xzs_ls, linterp, cmask

USE modd_prep_isba, ONLY : xgrid_soil, ngrid_level, lsnow_ideal,    &
                           xwsnow, xrsnow, xtsnow, xlwcsnow, xasnow,          &
                           xsg1snow, xsg2snow, xhistsnow, xagesnow


USE modd_isba_par,       ONLY : xwgmin
USE modd_data_cover_par, ONLY : nvegtype
USE modd_surf_par,       ONLY : xundef,nundef
!
USE modi_read_prep_isba_conf
USE modi_read_prep_isba_snow
USE modi_prep_isba_ascllv
USE modi_prep_isba_grib
USE modi_prep_isba_unif
USE modi_prep_isba_buffer
USE modi_abor1_sfx
USE modi_hor_interpol
USE modi_put_on_all_vegtypes
USE modi_vegtype_grid_to_patch_grid
USE modi_prep_hor_snow_fields
USE modi_get_luout
USE modi_prep_isba_extern
USE modi_prep_isba_netcdf
USE modi_vegtype_to_patch
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(data_cover_t), INTENT(INOUT) :: dtco
TYPE(isba_grid_t), INTENT(INOUT) :: ig
TYPE(isba_t), INTENT(INOUT) :: i
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  ! program calling surf. schemes
 CHARACTER(LEN=7),   INTENT(IN)  :: hsurf     ! type of field
 CHARACTER(LEN=28),  INTENT(IN)  :: hatmfile    ! name of the Atmospheric file
 CHARACTER(LEN=6),   INTENT(IN)  :: hatmfiletype! type of the Atmospheric file
 CHARACTER(LEN=28),  INTENT(IN)  :: hpgdfile    ! name of the Atmospheric file
 CHARACTER(LEN=6),   INTENT(IN)  :: hpgdfiletype! type of the Atmospheric file
!
LOGICAL, OPTIONAL,  INTENT(INOUT):: okey
!
!*      0.2    declarations of local variables
!
 CHARACTER(LEN=6)              :: yfiletype ! type of input file
 CHARACTER(LEN=28)             :: yfile     ! name of file
 CHARACTER(LEN=6)              :: yfiletype_snow ! type of input file
 CHARACTER(LEN=28)             :: yfile_snow     ! name of file
 CHARACTER(LEN=6)              :: yfilepgdtype_snow ! type of input file
 CHARACTER(LEN=28)             :: yfilepgd_snow     ! name of file 
 CHARACTER(LEN=6)              :: yfilepgdtype ! type of input file
 CHARACTER(LEN=28)             :: yfilepgd     ! name of file
REAL, POINTER, DIMENSION(:,:,:)     :: zfieldin  ! field to interpolate horizontally
REAL, POINTER, DIMENSION(:,:)       :: zfield, zpatch ! field to interpolate horizontally
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zfieldoutp ! field interpolated   horizontally
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zfieldoutv !
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zw        ! work array (x, fine   soil grid, npatch)
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zf        ! work array (x, output soil grid, npatch)
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zdg       ! out T grid (x, output soil grid, npatch)
INTEGER                       :: iluout    ! output listing logical unit
!
LOGICAL                       :: gunif     ! flag for prescribed uniform field
LOGICAL                       :: gunif_snow! flag for prescribed uniform field
INTEGER                       :: jpatch    ! loop on patches
INTEGER                       :: jvegtype  ! loop on vegtypes
INTEGER                       :: ini, inl, inp, jj, jl, ip_i, ip_o, jp, jveg ! Work integer
INTEGER, DIMENSION(SIZE(I%XDG,1),SIZE(I%XDG,3)) :: iwork
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------------
!
!
!*      1.     Reading of input file name and type
!
IF (lhook) CALL dr_hook('PREP_HOR_ISBA_FIELD',0,zhook_handle)
 CALL get_luout(hprogram,iluout)
!
 CALL read_prep_isba_conf(hprogram,hsurf,yfile,yfiletype,yfilepgd,yfilepgdtype,&
                         hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,iluout,gunif)
!
 cmask = 'NATURE'
!
ini=SIZE(ig%XLAT)
!
!-------------------------------------------------------------------------------------
!
!*      2.     Snow variables case?
!
IF (hsurf=='SN_VEG ') THEN
  CALL read_prep_isba_snow(hprogram,i%TSNOW%SCHEME,i%TSNOW%NLAYER,yfile_snow,yfiletype_snow,&
                                yfilepgd_snow,yfilepgdtype_snow,gunif_snow)
  IF(.NOT.gunif_snow.AND.len_trim(yfile_snow)==0.AND.len_trim(yfiletype_snow)==0)THEN
    IF(len_trim(yfile)/=0.AND.len_trim(yfiletype)/=0)THEN
       yfile_snow    =yfile
       yfiletype_snow=yfiletype
       yfilepgd_snow    =yfilepgd
       yfilepgdtype_snow=yfilepgdtype       
    ELSE
       gunif_snow=.true.
       IF(all(xwsnow==xundef))xwsnow=0.0
    ENDIF
  ENDIF
  CALL prep_hor_snow_fields(dtco, &
                            ig, u, &
                            hprogram, hsurf,                     &
                            yfile_snow, yfiletype_snow,          &
                            yfilepgd_snow, yfilepgdtype_snow,    &
                            iluout, gunif_snow, i%NPATCH, 1,     &
                            ini,i%TSNOW, i%TTIME,                &
                            xwsnow, xrsnow, xtsnow, xlwcsnow,    &
                            xasnow, lsnow_ideal, xsg1snow,       &
                            xsg2snow, xhistsnow, xagesnow,       &
                            i%XVEGTYPE, i%XVEGTYPE_PATCH, i%XPATCH,    &
                            okey                                 )
  DEALLOCATE(xwsnow)
  DEALLOCATE(xrsnow)
  DEALLOCATE(xtsnow)
  DEALLOCATE(xlwcsnow)
  DEALLOCATE(xsg1snow)
  DEALLOCATE(xsg2snow)
  DEALLOCATE(xhistsnow)
  DEALLOCATE(xagesnow)
  IF (lhook) CALL dr_hook('PREP_HOR_ISBA_FIELD',1,zhook_handle)
  RETURN
END IF
!
!-------------------------------------------------------------------------------------
!
!*      3.     Reading of input  configuration (Grid and interpolation type)
!
IF (gunif) THEN
  CALL prep_isba_unif(iluout,hsurf,zfieldin)
ELSE IF (yfiletype=='ASCLLV') THEN
  CALL prep_isba_ascllv(dtco, ug, u, uss, &
                        hprogram,hsurf,iluout,zfieldin)
ELSE IF (yfiletype=='GRIB  ') THEN
  CALL prep_isba_grib(hprogram,hsurf,yfile,iluout,zfieldin,okey)
ELSE IF (yfiletype=='MESONH' .OR. yfiletype=='ASCII ' .OR. yfiletype=='LFI   '.OR.yfiletype=='FA    ') THEN
   CALL prep_isba_extern(dtco, i, u, &
                         hprogram,hsurf,yfile,yfiletype,yfilepgd,yfilepgdtype,iluout,zfieldin,okey)
ELSE IF (yfiletype=='BUFFER') THEN
   CALL prep_isba_buffer(ig, u, &
                         hprogram,hsurf,iluout,zfieldin)
ELSE IF (yfiletype=='NETCDF') THEN
   CALL prep_isba_netcdf(dtco, u, &
                         hprogram,hsurf,yfile,iluout,zfieldin)
ELSE
   CALL abor1_sfx('PREP_HOR_ISBA_FIELD: data file type not supported : '//yfiletype)
END IF
!
!-------------------------------------------------------------------------------------
!
!*      5.     Horizontal interpolation
!
inl = SIZE(zfieldin,2)
inp = SIZE(zfieldin,3)
!
ALLOCATE(zfieldoutp(ini,inl,inp))
ALLOCATE(zfield(SIZE(zfieldin,1),inl))
!
ALLOCATE(zpatch(ini,inp))
zpatch(:,:) = 0.
!
IF (inp==nvegtype) THEN
  zpatch(:,:) = i%XVEGTYPE(:,:)
ELSEIF (inp==i%NPATCH) THEN
  zpatch(:,:) = i%XPATCH(:,:)
ELSEIF (inp<i%NPATCH) THEN
  DO jp = 1,i%NPATCH
    DO jveg = 1,nvegtype
      ip_i = vegtype_to_patch(jveg,inp)
      ip_o = vegtype_to_patch(jveg,i%NPATCH)
      !
      ! pour chaque patch d'entrée à interpoler, le masque
      ! est la somme des patchs de sortie (plus détaillés) présent sur 
      ! chaque point
      IF (ip_o==jp) THEN
        zpatch(:,ip_i) = zpatch(:,ip_i) + i%XPATCH(:,ip_o)
        EXIT
      ENDIF
    ENDDO
  ENDDO
ELSEIF (inp>i%NPATCH) THEN
  DO jp = 1,inp
    DO jveg = 1,nvegtype
      ip_i = vegtype_to_patch(jveg,inp)
      ip_o = vegtype_to_patch(jveg,i%NPATCH)
      !
      ! pour chaque patch d'entrée à interpoler, le masque
      ! est le patch de sortie (moins détaillé) présent 
      ! sur ce point
      IF (ip_i==jp) THEN
        zpatch(:,ip_i) = i%XPATCH(:,ip_o)
        EXIT
      ENDIF
    ENDDO
  ENDDO
ENDIF
!
DO jpatch = 1, inp

  zfield=zfieldin(:,:,jpatch)
  linterp(:) = (zpatch(:,jpatch) > 0.)
  CALL hor_interpol(dtco, u, &
                    iluout,zfield,zfieldoutp(:,:,jpatch))
  linterp = .true.

END DO
!
DEALLOCATE(zfield)
!
ALLOCATE(zfieldoutv(ini,inl,nvegtype))
!
 CALL put_on_all_vegtypes(ini,inl,inp,nvegtype,zfieldoutp,zfieldoutv)
!
DEALLOCATE(zfieldoutp)
!
!-------------------------------------------------------------------------------------
!
!*      6.     Transformation from vegtype grid to patch grid
!
ALLOCATE(zw(ini,SIZE(zfieldoutv,2),i%NPATCH))
!
zw = 0.
 CALL vegtype_grid_to_patch_grid(i%NPATCH,i%XVEGTYPE_PATCH,i%XPATCH,zfieldoutv,zw)
!
!-------------------------------------------------------------------------------------
!
!*      7.     Return to historical variable
!
!
SELECT CASE (hsurf)
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('ZS     ') 
  ALLOCATE(xzs_ls(ini))
  xzs_ls(:) = zfieldoutv(:,1,1)
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('WG     ') 
  ALLOCATE(zf(ini,i%NGROUND_LAYER,i%NPATCH))
  !
  !* interpolates on output levels
  CALL init_from_ref_grid(xgrid_soil,zw,i%XDG,zf)
  !
  !* retrieves soil water content from soil relative humidity
  ALLOCATE(i%XWG(ini,i%NGROUND_LAYER,i%NPATCH))
  i%XWG(:,:,:)=xundef
  IF(i%CISBA=='DIF')THEN
     iwork(:,:)=i%NWG_LAYER(:,:)
  ELSE
     iwork(:,:)=SIZE(i%XWG,2)
  ENDIF
  DO jpatch=1,i%NPATCH
    DO jj=1,ini
       IF(iwork(jj,jpatch)==nundef)cycle
       inl=iwork(jj,jpatch)
       DO jl=1,inl
          i%XWG(jj,jl,jpatch) = i%XWWILT(jj,jl) + zf(jj,jl,jpatch) * (i%XWFC(jj,jl)-i%XWWILT(jj,jl))
          i%XWG(jj,jl,jpatch) = max(min(i%XWG(jj,jl,jpatch),i%XWSAT(jj,jl)),xwgmin)
       ENDDO
    ENDDO
  ENDDO
  !
  WHERE(zf(:,:,:)==xundef)i%XWG(:,:,:)=xundef
  !
  DEALLOCATE(zf)
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('WGI    ')
  ALLOCATE(zf(ini,i%NGROUND_LAYER,i%NPATCH))
  !
  !* interpolates on output levels
  CALL init_from_ref_grid(xgrid_soil,zw,i%XDG,zf)
  !
  !* retrieves soil ice content from soil relative humidity
  ALLOCATE(i%XWGI(ini,i%NGROUND_LAYER,i%NPATCH))
  i%XWGI(:,:,:)=0.0
  IF(i%CISBA=='DIF')THEN
     iwork(:,:)=i%NWG_LAYER(:,:)
  ELSE
     iwork(:,:)=2
  ENDIF  
  DO jpatch=1,i%NPATCH
    DO jj=1,ini
       IF(iwork(jj,jpatch)==nundef)cycle
       inl=iwork(jj,jpatch)
       DO jl=1,inl
          i%XWGI(jj,jl,jpatch) = zf(jj,jl,jpatch) * i%XWSAT(jj,jl)
          i%XWGI(jj,jl,jpatch) = max(min(i%XWGI(jj,jl,jpatch),i%XWSAT(jj,jl)),0.)
       ENDDO
    ENDDO
  END DO
  !
  WHERE(zf(:,:,:)==xundef )i%XWGI(:,:,:)=xundef
  WHERE(i%XWGI(:,:,:)<=1.0e-10)i%XWGI(:,:,:)=0.0
  !
  DEALLOCATE(zf)
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('TG     ') 
  IF(i%LTEMP_ARP)THEN
    inl=i%NTEMPLAYER_ARP
  ELSE
    inl=i%NGROUND_LAYER
  ENDIF
  ALLOCATE(i%XTG(ini,inl,i%NPATCH))
  ALLOCATE(zdg(SIZE(i%XDG,1),inl,SIZE(i%XDG,3)))
  IF (i%CISBA=='2-L'.OR.i%CISBA=='3-L') THEN
     DO jpatch=1,i%NPATCH
        zdg(:,1,jpatch) = 0.01
        zdg(:,2,jpatch) = 0.40                    ! deep temperature for force-restore taken at 20cm
        IF(i%CISBA=='3-L') zdg(:,3,jpatch) = 5.00   ! climatological temperature, usually not used
     ENDDO         
     IF(i%LTEMP_ARP)THEN
       DO jpatch=1,i%NPATCH
          zdg(:,3,jpatch) = 1.0
          DO jl=4,inl
             zdg(:,jl,jpatch) = zdg(:,jl-1,jpatch)+1.0
          ENDDO
       ENDDO
     ENDIF
  ELSE
    !* diffusion method, the soil grid is the same as for humidity
    zdg(:,:,:) = i%XDG(:,:,:)
  END IF
  CALL init_from_ref_grid(xgrid_soil,zw,zdg,i%XTG)
  DEALLOCATE(zdg)
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('WR     ') 
  ALLOCATE(i%XWR(ini,i%NPATCH))
  DO jpatch=1,i%NPATCH
    i%XWR(:,jpatch) = zw(:,1,jpatch)
  END DO
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
  CASE('WRL    ') 
  ALLOCATE(i%XWRL(ini,i%NPATCH))
  DO jpatch=1,i%NPATCH
    i%XWRL(:,jpatch) = zw(:,1,jpatch)
  END DO
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
  CASE('WRLI   ') 
  ALLOCATE(i%XWRLI(ini,i%NPATCH))
  DO jpatch=1,i%NPATCH
    i%XWRLI(:,jpatch) = zw(:,1,jpatch)
  END DO
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('WRVN   ') 
  ALLOCATE(i%XWRVN(ini,i%NPATCH))
  DO jpatch=1,i%NPATCH
    i%XWRVN(:,jpatch) = zw(:,1,jpatch)
  END DO
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('TV     ') 
  ALLOCATE(i%XTV(ini,i%NPATCH))
  DO jpatch=1,i%NPATCH
    i%XTV(:,jpatch) = zw(:,1,jpatch)
  END DO
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('TL     ') 
  ALLOCATE(i%XTL(ini,i%NPATCH))
  DO jpatch=1,i%NPATCH
    i%XTL(:,jpatch) = zw(:,1,jpatch)
  END DO
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('TC     ') 
  ALLOCATE(i%XTC(ini,i%NPATCH))
  DO jpatch=1,i%NPATCH
    i%XTC(:,jpatch) = zw(:,1,jpatch)
  END DO
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('QC     ') 
  ALLOCATE(i%XQC(ini,i%NPATCH))
  DO jpatch=1,i%NPATCH
    i%XQC(:,jpatch) = zw(:,1,jpatch)
  END DO
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('LAI    ') 
  !* LAI is updated only if present and pertinent (evolutive LAI) in input file
  IF (any(zw(:,:,:)/=xundef)) THEN
    DO jpatch=1,i%NPATCH
      i%XLAI(:,jpatch) = zw(:,1,jpatch)
    END DO
  END IF
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('ICE_STO') 
  ALLOCATE(i%XICE_STO(ini,i%NPATCH))
  DO jpatch=1,i%NPATCH
    i%XICE_STO(:,jpatch) = zw(:,1,jpatch)
  END DO
  !
END SELECT
!
DEALLOCATE(zw)
!-------------------------------------------------------------------------------------
!
!*      8.     Deallocations
!
DEALLOCATE(zfieldin )
DEALLOCATE(zfieldoutv)
!
IF (lhook) CALL dr_hook('PREP_HOR_ISBA_FIELD',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------
!
 CONTAINS
!
!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------
!
SUBROUTINE init_from_ref_grid(PGRID1,PT1,PD2,PT2)
!
USE modi_interp_grid_nat
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pt1    ! variable profile
REAL, DIMENSION(:),     INTENT(IN)  :: pgrid1 ! normalized grid
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pd2    ! output layer thickness
REAL, DIMENSION(:,:,:), INTENT(OUT) :: pt2    ! variable profile
!
INTEGER                                  :: ji,jl  ! loop counter
REAL, DIMENSION(SIZE(PT1,1),SIZE(PT1,2)) :: zd1 ! input grid
REAL, DIMENSION(SIZE(PD2,1),SIZE(PD2,2)) :: zd2 ! output grid
!
INTEGER :: ilayer1, ilayer2
REAL(KIND=JPRB) :: zhook_handle
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (lhook) CALL dr_hook('INIT_FROM_REF_GRID',0,zhook_handle)
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IF (SIZE(pt1,2)==3) THEN
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!* 1. case with only 3 input levels (typically coming from 'UNIF')
!     -----------------------------
!
  IF (i%CISBA=='2-L' .OR. i%CISBA=='3-L') THEN
    !* Possible LTEMP_ARP case
    IF(SIZE(pt2,2)>3)THEN
      ilayer1=3
      ilayer2=SIZE(pt2,2)
    ELSE
      ilayer1=SIZE(pt2,2)
      ilayer2=0
    ENDIF
    !* historical 2L or 3L ISBA version
    DO jpatch=1,i%NPATCH
      pt2(:,1:ilayer1,jpatch) = pt1(:,1:ilayer1,jpatch) 
      !* Possible LTEMP_ARP case
      IF(ilayer2>0)THEN
        DO jl=ilayer1+1,ilayer2
          pt2(:,jl,jpatch) = pt2(:,ilayer1,jpatch)
        ENDDO
      ENDIF
    END DO
  ELSEIF(i%CISBA=='DIF')THEN
    DO jpatch=1,i%NPATCH
       !surface layer (generally 0.01m imposed)
       pt2(:,1,jpatch) = pt1(:,1,jpatch) 
       !second layer
       pt2(:,2,jpatch) = 0.25*pt1(:,1,jpatch)+0.75*pt1(:,2,jpatch)
       !others layers
       DO ji=1,SIZE(pt1,1)
         DO jl=3,i%NGROUND_LAYER
           IF(pd2(ji,jl,jpatch)<=i%XDG2(ji,jpatch))THEN 
            !root layers
             pt2(ji,jl,jpatch) = pt1(ji,2,jpatch)
           ELSE
            !deep layers
             pt2(ji,jl,jpatch) = pt1(ji,3,jpatch)
           ENDIF
         END DO
       END DO 
    END DO 
  END IF    
!  
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
ELSE
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!* 2. case with fine grid as input (general case)
!     ----------------------------
!
  DO jl=1,SIZE(pt1,2)
     zd1(:,jl) = pgrid1(jl)
  ENDDO
!
  DO jpatch=1,i%NPATCH
     zd2(:,:) = pd2(:,:,jpatch)
     CALL interp_grid_nat(zd1,pt1(:,:,jpatch),zd2,pt2(:,:,jpatch))
  ENDDO
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
ENDIF
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (lhook) CALL dr_hook('INIT_FROM_REF_GRID',1,zhook_handle)
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
END SUBROUTINE init_from_ref_grid
!-------------------------------------------------------------------------------------
!
END SUBROUTINE prep_hor_isba_field