prep_trip_run.F90

Go to the documentation of this file.

!     #########
      SUBROUTINE prep_trip_run (TP, TPG, &
                                KYEAR,KMONTH,KDAY,PTIME,KLON,KLAT)
!     ####################
!
!!****  *PREP_TRIP_RUN*  
!!
!!    PURPOSE
!!    -------
!
!     Prepare TRIP variables and parameters.
!     
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!      
!!    REFERENCE
!!    ---------
!!      
!!    AUTHOR
!!    ------
!!      B. Decharme     (Meteo-France)
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    27/05/08 
!!      09/16   B. Decharme  limit wtd to -1000m
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
USE modd_trip, ONLY : trip_t
USE modd_trip_grid, ONLY : trip_grid_t
!
USE modn_trip, ONLY : cgroundw, cvit, lflood,  &
                      xcvel, xratmed, xtstep,  &
                      lgwsubf, xgwsubd

USE modn_trip_prep, ONLY : xtaug_unif, lgweq, lread_flood, &
                           xtaug_up, xtaug_down
!
USE modd_trip_par
USE modd_trip_listing, ONLY : nlisting
!
USE mode_trip_init
USE mode_rw_trip
!
USE modi_abort_trip
USE modi_get_trip_grid_conf
USE modi_init_param_trip
USE modi_init_restart_trip
USE modi_get_lonlat_trip
!
USE modi_gwf
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!-------------------------------------------------------------------------------
!
!*      0.1    declarations of arguments
!
!
TYPE(trip_t), INTENT(INOUT) :: TP
TYPE(trip_grid_t), INTENT(INOUT) :: TPG
!
INTEGER,          INTENT(IN) :: KYEAR   !date UTC
INTEGER,          INTENT(IN) :: KMONTH  !date UTC
INTEGER,          INTENT(IN) :: KDAY    !date UTC
REAL,             INTENT(IN) :: PTIME   !date UTC
INTEGER,          INTENT(OUT):: KLON    ! number of points in longitude
INTEGER,          INTENT(OUT):: KLAT    ! number of points in latitude
!
!-------------------------------------------------------------------------------
!
!*      0.2    declarations of local variables
!
LOGICAL,           PARAMETER         :: LDOUBLE         =.true.
!
REAL,              PARAMETER         :: ZMINCELL        = 3.0
!
 CHARACTER(LEN=20), PARAMETER         :: YFILE_PARAM_1D  ='TRIP_PGD_1D.nc'
 CHARACTER(LEN=20), PARAMETER         :: YFILE_PARAM_HD  ='TRIP_PGD_HD.nc'
 CHARACTER(LEN=20), PARAMETER         :: YFILE_PARAM_12  ='TRIP_PGD_12D.nc'
 CHARACTER(LEN=20), PARAMETER         :: YFILE_INIT_1D   ='TRIP_INIT_1D.nc'
 CHARACTER(LEN=20), PARAMETER         :: YFILE_INIT_HD   ='TRIP_INIT_HD.nc'
 CHARACTER(LEN=20), PARAMETER         :: YFILE_INIT_12   ='TRIP_INIT_12D.nc'
!
 CHARACTER(LEN=20), PARAMETER         :: YFILE_PARAM     ='TRIP_PARAM.nc'
 CHARACTER(LEN=20), PARAMETER         :: YFILE_PREP      ='TRIP_PREP.nc'
!
 CHARACTER(LEN=20)                    :: YFILE_READ
 CHARACTER(LEN=20)                    :: YFILE_READ_INIT
!
 CHARACTER(LEN=20), PARAMETER         :: YFILE_FLOOD_1D  ='TRIP_SGFLOOD_1D.nc'
 CHARACTER(LEN=20), PARAMETER         :: YFILE_FLOOD_HD  ='TRIP_SGFLOOD_HD.nc'
 CHARACTER(LEN=20)                    :: YFILE_FLOOD_READ
!
 CHARACTER(LEN=20), PARAMETER         :: YFILE_GW_HD     ='TRIP_SGAQUI_HD.nc'
 CHARACTER(LEN=20), PARAMETER         :: YFILE_GW_12D    ='TRIP_SGAQUI_12D.nc'
 CHARACTER(LEN=20)                    :: YFILE_GW_READ
 CHARACTER(LEN=20), PARAMETER         :: YFILE_GW_FRC_HD ='GWEQ_HD_FRC.nc'
!
 CHARACTER(LEN=20)                    :: YVAR 
!
REAL, DIMENSION(:),ALLOCATABLE       :: ZLON
REAL, DIMENSION(:),ALLOCATABLE       :: ZLAT
!
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZGRCN
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZSEQ
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZFRACAREA
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZNUM_BAS
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZDR_AREA
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZREAD
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZWORK
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZLEN
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZGREEN_ANT
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZHSTREAM
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZDRAIN
REAL,DIMENSION(:,:),ALLOCATABLE      :: ZHG_OLD
!
REAL,DIMENSION(:,:,:),ALLOCATABLE    :: ZREAD3D
!
INTEGER,DIMENSION(:,:),ALLOCATABLE   :: IGRCN
INTEGER,DIMENSION(:,:),ALLOCATABLE   :: INEXTX
INTEGER,DIMENSION(:,:),ALLOCATABLE   :: INEXTY
!
INTEGER,DIMENSION(:), ALLOCATABLE    :: INCELL_BAS
!
REAL    :: ZLONMIN   ! minimum longitude (degrees)
REAL    :: ZLONMAX   ! maximum longitude (degrees)
REAL    :: ZLATMIN   ! minimum latitude  (degrees)
REAL    :: ZLATMAX   ! maximum latitude  (degrees)
REAL    :: ZGRID_RES ! 1° or 0.5° resolution
REAL    :: ZSEQMAX
!
INTEGER :: ILON,ILAT
INTEGER :: ILON_G, ILON_DEB, ILON_END
INTEGER :: ILON_G12, ILON_DEB12, ILON_END12
INTEGER :: ILAT_G, ILAT_DEB, ILAT_END
INTEGER :: ILAT_G12, ILAT_DEB12, ILAT_END12
!
INTEGER :: IWORK, IFLOOD, INI, ILATF
INTEGER :: JLON, JLAT, JBAS, JITER
!
LOGICAL :: LMASKLON, LMASKLAT, LWORK
!
!-------------------------------------------------------------------------------
!Output attribut for netcdf diag file
!-------------------------------------------------------------------------------
!
 CHARACTER(LEN=50) :: YTITLE,YTIMEUNIT
!                    YTITLE    = Title of each output file
!                    YTIMEUNIT = Time unit in each output file if present
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!Initilyse TRIP
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('PREP_TRIP_RUN',0,zhook_handle)
!
WRITE(nlisting,*)''
WRITE(nlisting,*)'        PREP TRIP            '
WRITE(nlisting,*)''
!
!-------------------------------------------------------------------------------
!
IF(lflood)THEN
  IF(cgroundw=='DEF') &
      WRITE(nlisting,*)'! Attention, you use the flooding scheme without groundwater scheme !!!'       
  IF(cvit /= 'VAR')THEN
    WRITE(nlisting,*)'! You cannot use the flooding scheme without the variable velocity scheme !!!'
    CALL abort_trip('PREP_TRIP_RUN: You cannot use the flooding scheme without the variable velocity scheme')
  ENDIF
ENDIF
!
IF(cgroundw=='CST')THEN
  IF(xtaug_unif<1.0.OR.xtaug_unif>365.0)THEN
    WRITE(nlisting,*)'! Constant transfert time value XTAUG_UNIF must be at least 1 day or inferior to 365 days !!!'
    CALL abort_trip('PREP_TRIP_RUN: Constant transfert time value must be at least 1 day or inferior to 365 days')
  ENDIF        
ENDIF
!
IF(cgroundw=='DIF')THEN
  IF(cvit /= 'VAR')THEN
    WRITE(nlisting,*)'! You cannot use the groundwater scheme without the variable velocity scheme !!!'
    CALL abort_trip('PREP_TRIP_RUN: You cannot use the groundwater scheme without the variable velocity scheme')
  ENDIF
  IF(xtaug_up<1.0)THEN
    WRITE(nlisting,*)'! Upstream transfert time value XTAUG_UP must be at least 1 day !!!'
    CALL abort_trip('PREP_TRIP_RUN: Upstream transfert time value must be at least 1 day')
  ENDIF
  IF(xtaug_down>365.0)THEN
    WRITE(nlisting,*)'! Downstream transfert time value XTAUG_DOWN must be lower than 365 days !!!'
    CALL abort_trip('PREP_TRIP_RUN: Downstream transfert time value must be lower than 365 days')
  ENDIF  
ENDIF
!
IF(lgwsubf.AND.xgwsubd>30.)THEN
     WRITE(nlisting,*)'!! XGWSUBD too large (must be <=30), check your namelist  !!!' 
     CALL abort_trip('PREP_TRIP_RUN: XGWSUBD too large (must be <=30), check your namelist !!!')
ENDIF
!
!-------------------------------------------------------------------------------
! * Read TRIP grid configuration
!-------------------------------------------------------------------------------
!
 CALL get_trip_grid_conf(tpg, &
                        zlonmin,zlonmax,zlatmin,zlatmax,zgrid_res,ilon,ilat)
klon=ilon
klat=ilat
!
ALLOCATE(zlon(ilon))
ALLOCATE(zlat(ilat))
zlon(:)=xundef
zlat(:)=xundef
!
 CALL get_lonlat_trip(tpg, &
                     klon,klat,zlon,zlat)
!
IF(zgrid_res==1.0)THEN
!
  WRITE(nlisting,*)'! 1° by 1° TRIP run !!!'
!
  ilon_g = 360
  ilat_g = 180
  ilatf  = 31
  yfile_read = yfile_param_1d
  yfile_read_init = yfile_init_1d
  yfile_flood_read = yfile_flood_1d
!
ELSEIF(zgrid_res==0.5)THEN
!
  WRITE(nlisting,*)'! 0.5° by 0.5° TRIP run !!!'
!
  ilon_g = 720
  ilat_g = 360
  ilatf  = 61
  yfile_read = yfile_param_hd
  yfile_read_init = yfile_init_hd
  yfile_gw_read = yfile_gw_hd
  yfile_flood_read = yfile_flood_hd
!
ELSEIF(zgrid_res < 0.09)THEN
!
  WRITE(nlisting,*)'!!! 1/12° by 1/12° TRIP run (FRANCE) !!!'
!
  ilon_g = 180
  ilat_g = 108
  ilon_g12 = 4320
  ilat_g12 = 2160
  
  ilatf  = 1
  yfile_read = yfile_param_12
  yfile_read_init = yfile_init_12
  yfile_gw_read = yfile_gw_12d  
!
  IF(xratmed==1.4)THEN
     WRITE(nlisting,*)'! meandering ratio is the same than at 1° resolution   !!!'   
     WRITE(nlisting,*)'! change XRATMED in namelist                           !!!'           
     CALL abort_trip('PREP_TRIP_RUN: meandering ratio not good')
  ENDIF 
!
ELSE
!
  WRITE(nlisting,*)'! The resolution of the TRIP run is not good !!!'   
  WRITE(nlisting,*)'! Should be 1° or 0.5° or 1/12 ° over France !!!'
  CALL abort_trip('PREP_TRIP_RUN: resolution of the TRIP run is not good')
!  
ENDIF
!
IF(cgroundw=='DIF'.AND.zgrid_res>0.5)THEN
  WRITE(nlisting,*)'! You cannot use the groundwater scheme with another resolution than 0.5 or 1/12 !!!'
  CALL abort_trip('PREP_TRIP_RUN: You cannot use the groundwater scheme with another resolution than 0.5 or 1/12')
ENDIF  
!
IF(lgweq.AND.zgrid_res/=0.5)THEN
  WRITE(nlisting,*)'! You cannot use groundwater equilibrium with another resolution than 0.5  !!!'
  CALL abort_trip('PREP_TRIP_RUN: You cannot use groundwater equilibrium with another resolution than 0.5')
ENDIF
!
lmasklon=.false.
IF(zlonmin/=-180..OR.zlonmax/=180.)lmasklon=.true.
!
lmasklat=.false.
IF(zlatmin/=-180..OR.zlatmax/=180.)lmasklat=.true.
!
!-------------------------------------------------------------------------------
!Allocate arguments
!-------------------------------------------------------------------------------
!
ALLOCATE(tpg%GMASK    (ilon,ilat))
ALLOCATE(tpg%GMASK_VEL(ilon,ilat))
ALLOCATE(tpg%GMASK_GW (ilon,ilat))
ALLOCATE(tpg%GMASK_FLD(ilon,ilat))
ALLOCATE(tpg%GMASK_GRE(ilon,ilat))
ALLOCATE(tpg%GMASK_ANT(ilon,ilat))
tpg%GMASK    (:,:) = .false.
tpg%GMASK_VEL(:,:) = .false.
tpg%GMASK_GW (:,:) = .false.
tpg%GMASK_FLD(:,:) = .false.
tpg%GMASK_GRE(:,:) = .false.
tpg%GMASK_ANT(:,:) = .false.
!
ALLOCATE(tpg%XAREA (ilon,ilat))
ALLOCATE(tpg%NGRCN (ilon,ilat))
ALLOCATE(tpg%NBASID(ilon,ilat))
!
ALLOCATE(zlen(ilon,ilat))
ALLOCATE(zgrcn(ilon,ilat))       
ALLOCATE(zseq(ilon,ilat))   
ALLOCATE(zfracarea(ilon,ilat))
ALLOCATE(znum_bas(ilon,ilat))
ALLOCATE(zdr_area(ilon,ilat))
ALLOCATE(zgreen_ant(ilon,ilat))
!
ALLOCATE(tp%XSURF_STO(ilon,ilat))
!
IF(cgroundw/='DEF')THEN
  ALLOCATE(tp%XTAUG      (ilon,ilat))
ELSE
  ALLOCATE(tp%XTAUG      (0,0))
ENDIF
!
IF(cgroundw=='CST')THEN
  ALLOCATE(tp%XGROUND_STO(ilon,ilat))      
ELSE
  ALLOCATE(tp%XGROUND_STO(0,0))      
ENDIF
!
IF(cgroundw=='DIF')THEN
  ALLOCATE(tp%XHGROUND   (ilon,ilat))
  ALLOCATE(tp%XWEFF      (ilon,ilat))
  ALLOCATE(tp%XTRANS     (ilon,ilat))
  ALLOCATE(tp%XNUM_AQUI  (ilon,ilat))
  ALLOCATE(tp%XTOPO_RIV  (ilon,ilat))
ELSE
  ALLOCATE(tp%XHGROUND   (0,0))
  ALLOCATE(tp%XWEFF      (0,0))
  ALLOCATE(tp%XTRANS     (0,0))
  ALLOCATE(tp%XNUM_AQUI  (0,0))
  ALLOCATE(tp%XTOPO_RIV  (0,0))
ENDIF
!
IF(lflood.OR.cgroundw=='DIF')THEN
  ALLOCATE(tp%XHC_BED    (ilon,ilat))
ELSE
  ALLOCATE(tp%XHC_BED   (0,0))
ENDIF
!
IF(lflood)THEN
  ALLOCATE(tp%XN_FLOOD    (ilon,ilat))
  ALLOCATE(tp%XFLOOD_STO  (ilon,ilat))
  ALLOCATE(tp%XHFLOOD     (ilon,ilat))
  ALLOCATE(tp%XFFLOOD     (ilon,ilat))
ELSE
  ALLOCATE(tp%XN_FLOOD  (0,0))
  ALLOCATE(tp%XFLOOD_STO(0,0))
  ALLOCATE(tp%XHFLOOD   (0,0))
  ALLOCATE(tp%XFFLOOD   (0,0))
ENDIF
!
IF(cvit=='VAR')THEN
  ALLOCATE(tp%XSLOPEBED(ilon,ilat))      
  ALLOCATE(tp%XWIDTH   (ilon,ilat))      
  ALLOCATE(tp%XN       (ilon,ilat))
ELSE
  ALLOCATE(tp%XSLOPEBED(0,0))      
  ALLOCATE(tp%XWIDTH   (0,0))      
  ALLOCATE(tp%XN       (0,0))
ENDIF
!
!-------------------------------------------------------------------------------
! * Create param and init file
!-------------------------------------------------------------------------------
!
ytitle    = 'TRIP parameters for a run'
ytimeunit = '-'
 CALL init_param_trip(tpg, &
                     nlisting,yfile_param,ilon,ilat,ytitle,ytimeunit)
!
ytitle   ='TRIP prep historical variable'
ytimeunit='-'
 CALL init_restart_trip(tpg, &
                       nlisting,yfile_prep,ilon,ilat,ytitle,ytimeunit,.false.)
!
!-------------------------------------------------------------------------------
! * Compute the mask of the run
!-------------------------------------------------------------------------------        
!
IF(zgrid_res > 0.09)THEN      
!        
  ilon_deb = int(zlonmin/zgrid_res) + ilon_g/2 + 1
  ilon_end = int(zlonmax/zgrid_res) + ilon_g/2  
  ilat_deb = int(zlatmin/zgrid_res) + ilat_g/2 + 1
  ilat_end = int(zlatmax/zgrid_res) + ilat_g/2
!  
  ALLOCATE(igrcn(ilon_g,ilat_g))
  ALLOCATE(inextx(ilon_g,ilat_g))
  ALLOCATE(inexty(ilon_g,ilat_g))
  ALLOCATE(zwork(ilon_g,ilat_g))
!  
ELSE
!        
  ilon_deb = 1
  ilon_end = ilon_g
  ilat_deb = 1
  ilat_end = ilat_g
!        
  ilon_deb12 = int(zlonmin/zgrid_res) + ilon_g12/2 + 1
  ilon_end12 = int(zlonmax/zgrid_res) + ilon_g12/2  
  ilat_deb12 = int(zlatmin/zgrid_res) + ilat_g12/2 + 1
  ilat_end12 = int(zlatmax/zgrid_res) + ilat_g12/2
!  
  ALLOCATE(igrcn(ilon_g12,ilat_g12))
  ALLOCATE(inextx(ilon_g12,ilat_g12))
  ALLOCATE(inexty(ilon_g12,ilat_g12))
  ALLOCATE(zwork(ilon_g12,ilat_g12))
!  
ENDIF
!
ALLOCATE(zread(ilon_g,ilat_g))
!
!-------------------------------------------------------------------------------
! * Compute and Read TRIP parameter
!-------------------------------------------------------------------------------        
!
igrcn =0
inextx=0
inexty=0
!
! * Basin ID
!
yvar ='NUM_BAS'
 CALL read_trip(nlisting,yfile_read,yvar,zread)
znum_bas(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
!
WHERE(znum_bas>=xundef-1.0)znum_bas=0.0
tpg%NBASID(:,:)=int(znum_bas(:,:))
!
tpg%NBASMIN = minval(tpg%NBASID(:,:),tpg%NBASID(:,:)>0)
tpg%NBASMAX = maxval(tpg%NBASID(:,:),tpg%NBASID(:,:)>0)
!
! * Drainage Area
!
yvar ='DR_AREA'
 CALL read_trip(nlisting,yfile_read,yvar,zread)
zdr_area(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
!
! * Flow direction
!
yvar ='FLOWDIR'
 CALL read_trip(nlisting,yfile_read,yvar,zread)
zgrcn(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
WHERE(zread(:,:)==xundef)zread(:,:)=0.0
!
IF(zgrid_res > 0.09)THEN      
!
  igrcn = int(zread)
!
! * Set the distance between grids with the meandering ratio
!
  CALL setnext(ilon_g,ilat_g,igrcn,inextx,inexty)
!
  CALL setlen(ilon_g,ilat_g,igrcn,inextx,inexty,xratmed,zwork)
  zlen(:,:) = zwork(ilon_deb:ilon_end,ilat_deb:ilat_end)
!
ELSE
!
  igrcn(ilon_deb12:ilon_end12,ilat_deb12:ilat_end12) = int(zread)
!
! * Set the distance between grids with the meandering ratio
!
  CALL setnext(ilon_g12,ilat_g12,igrcn,inextx,inexty)
!
  CALL setlen(ilon_g12,ilat_g12,igrcn,inextx,inexty,xratmed,zwork)
  zlen(:,:) = zwork(ilon_deb12:ilon_end12,ilat_deb12:ilat_end12)
!
ENDIF
!
!
! * Set area size
!
IF(zgrid_res==0.5)THEN
  yvar ='FRACAREA'
  CALL read_trip(nlisting,yfile_read,yvar,zread)
  zfracarea(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
ELSE
  zfracarea(:,:) = 1.0
ENDIF
!
 CALL setarea(ilat,zlatmin,zgrid_res,tpg%XAREA)
!
zlen(:,:) = zlen(:,:) * zfracarea(:,:)
tpg%XAREA(:,:) = tpg%XAREA(:,:) * zfracarea(:,:)
!
DEALLOCATE(igrcn )
DEALLOCATE(inextx)
DEALLOCATE(inexty)
DEALLOCATE(zwork )
DEALLOCATE(zfracarea)   
!
! * Domain test
!
ALLOCATE(inextx(ilon,ilat))
ALLOCATE(inexty(ilon,ilat))
!
WHERE(zgrcn(:,:)<xundef-1.0)
      tpg%NGRCN(:,:)=int(zgrcn(:,:))
      tpg%GMASK(:,:)=.true.
ELSEWHERE
      tpg%NGRCN(:,:)=0
      tpg%GMASK(:,:)=.false.
ENDWHERE
!
 CALL setnext(ilon,ilat,tpg%NGRCN,inextx,inexty,lmasklon,lmasklat)
!
! * Store some param
!
yvar ='FLOWDIR'
zgrcn(:,:)=REAL(tpg%ngrcn(:,:))
 CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,zgrcn)
yvar ='RIVLEN'
 CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,zlen,odouble=ldouble)
yvar ='NUM_BAS'
 CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,znum_bas)
yvar ='DR_AREA'
 CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,zdr_area,odouble=ldouble)
yvar ='CELL_AREA'
 CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,tpg%XAREA,odouble=ldouble)
!
DEALLOCATE(inextx)
DEALLOCATE(inexty)
DEALLOCATE(zdr_area)
DEALLOCATE(znum_bas)
!
! * River sequence values read
!
yvar ='RIVSEQ'
 CALL read_trip(nlisting,yfile_read,yvar,zread)
zseq(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
!
 CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,zseq)
!
! * Initial river storage
!
yvar ='SURF_STO'
 CALL read_trip(nlisting,yfile_read_init,yvar,zread)
tp%XSURF_STO(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
!        
 CALL write_trip(nlisting,yfile_prep,yvar,tpg%GMASK,tp%XSURF_STO,odouble=ldouble)
!
! * Greenland and Antarctic masks
!
yvar ='GREEN_ANT'
IF(zgrid_res > 0.09)THEN
  CALL read_trip(nlisting,yfile_read,yvar,zread)
ELSE
  zread(:,:)=0.0
ENDIF
zgreen_ant(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
 CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,zgreen_ant)
DO jlat=1,ilat
   DO jlon=1,ilon
      IF(zgreen_ant(jlon,jlat)==2.0)THEN
          tpg%GMASK_ANT(jlon,jlat)=.true.
      ELSEIF(zgreen_ant(jlon,jlat)==1.0)THEN
          tpg%GMASK_GRE(jlon,jlat)=.true.
      ENDIF
   ENDDO
ENDDO
!
DEALLOCATE(zlon)
DEALLOCATE(zlat)
DEALLOCATE(zgreen_ant)
!
! * Variable velocity scheme parameters
!
IF(cvit == 'VAR')THEN
!
  yvar ='N_RIV'
  CALL read_trip(nlisting,yfile_read,yvar,zread)
  tp%XN(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,tp%XN)
!
  yvar ='WIDTHRIV'
  CALL read_trip(nlisting,yfile_read,yvar,zread)
  tp%XWIDTH(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,tp%XWIDTH)
!
  yvar ='SLOPERIV'
  CALL read_trip(nlisting,yfile_read,yvar,zread)
  tp%XSLOPEBED(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
  WHERE(tp%XSLOPEBED(:,:)<xundef-1.0)
        tp%XSLOPEBED(:,:) = max(1.e-5,tp%XSLOPEBED(:,:)/xratmed)
  ENDWHERE
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,tp%XSLOPEBED,odouble=ldouble)  
!  
  tpg%GMASK_VEL(:,:)=tpg%GMASK(:,:)
  WHERE(tp%XWIDTH(:,:)>=xundef-1.0)
        tpg%GMASK_VEL(:,:)=.false.
  ENDWHERE
!  
ENDIF 
!
! * River depth
!
IF (lflood.OR.cgroundw == 'DIF') THEN
!
  yvar ='RIVDEPTH'     
  CALL read_trip(nlisting,yfile_read,yvar,zread)
  tp%XHC_BED(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
  WHERE(tp%XHC_BED(:,:)==0.0)
        tp%XHC_BED(:,:)=xundef
  ENDWHERE
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK,tp%XHC_BED)
!
ENDIF  
!
! * Groundwater parameters
!
IF(cgroundw=='CST')THEN
!
  tpg%GMASK_GW(:,:)=tpg%GMASK(:,:)
!
  yvar ='GROUND_STO'
  CALL read_trip(nlisting,yfile_read_init,yvar,zread)
  tp%XGROUND_STO(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
!  
  tp%XTAUG(:,:)=xtaug_unif
  WHERE(tpg%GMASK_ANT(:,:).OR.tpg%GMASK_GRE(:,:))
        tp%XTAUG      (:,:)=xundef
        tpg%GMASK_GW  (:,:)=.false.
        tp%XGROUND_STO(:,:)=0.0
  ENDWHERE
!          
  yvar ='TAUG'
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_GW,tp%XTAUG)
!  
  yvar ='GROUND_STO'
  CALL write_trip(nlisting,yfile_prep,yvar,tpg%GMASK_GW,tp%XGROUND_STO,odouble=ldouble)
!
ELSEIF(cgroundw=='DIF')THEN
!
  tpg%GMASK_GW(:,:)=tpg%GMASK_VEL(:,:)
!
  yvar ='NUM_AQUI'
  CALL read_trip(nlisting,yfile_read,yvar,zread)
  tp%XNUM_AQUI(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
!  
  WHERE(tp%XNUM_AQUI(:,:)>=xundef-1.0)
       tpg%GMASK_GW(:,:)=.false.
  ENDWHERE
!  
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_GW,tp%XNUM_AQUI)  
!
  IF(zgrid_res > 0.09)THEN
    zseqmax=5.0
  ELSE
    zseqmax=10.0
  ENDIF
  tp%XTAUG(:,:)=xtaug_down+(xtaug_up-xtaug_down)*(zseqmax-min(zseqmax,zseq(:,:)))/(zseqmax-1.0)
!  
  yvar ='TAUG'
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_GW,tp%XTAUG)
!
  yvar ='WEFF'
  CALL read_trip(nlisting,yfile_read,yvar,zread)
  tp%XWEFF(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_GW,tp%XWEFF)
!
  yvar ='TRANS'
  CALL read_trip(nlisting,yfile_read,yvar,zread)
  tp%XTRANS(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_GW,tp%XTRANS)
!
  yvar ='TOPO_RIV'  
  CALL read_trip(nlisting,yfile_read,yvar,zread)  
  tp%XTOPO_RIV(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_GW,tp%XTOPO_RIV) 
!
  ALLOCATE(zread3d(ilon_g,ilat_g,ndimtab))
  zread3d(:,:,:)= xundef
!
  ALLOCATE(tp%XTABGW_F (ilon,ilat,ndimtab))      
  ALLOCATE(tp%XTABGW_H (ilon,ilat,ndimtab))      
!  
  tp%XTABGW_F(:,:,:)= xundef
  tp%XTABGW_H(:,:,:)= xundef
!
  yvar ='TABGW_F'
  CALL read_trip(nlisting,yfile_gw_read,yvar,zread3d(:,ilatf:ilat_g,:))
  tp%XTABGW_F(:,:,:) = zread3d(ilon_deb:ilon_end,ilat_deb:ilat_end,:)
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_GW,tp%XTABGW_F)
!
  yvar ='TABGW_H'
  CALL read_trip(nlisting,yfile_gw_read,yvar,zread3d(:,ilatf:ilat_g,:))
  tp%XTABGW_H(:,:,:) = zread3d(ilon_deb:ilon_end,ilat_deb:ilat_end,:)
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_GW,tp%XTABGW_H)
!
  DEALLOCATE(zread3d)
!
! Groundwater init from file
!          
  yvar ='HGROUND'   
  CALL read_trip(nlisting,yfile_read_init,yvar,zread)   
  tp%XHGROUND(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end) 
!
  WHERE(tpg%GMASK_GW(:,:))
        tp%XHGROUND(:,:) = max(tp%XHGROUND(:,:),tp%XTOPO_RIV(:,:)-xgwdzmax)
  ENDWHERE
!
  IF(lgweq)THEN
!
!   Groundwater equilibrium (appendix B in Vergnes thesis)
!
    ALLOCATE(zhstream(ilon,ilat))   
    ALLOCATE(zdrain(ilon,ilat))  
    ALLOCATE(zhg_old(ilon,ilat))
!
    lwork         = .false.
    zhg_old(:,:) = xundef
!
    yvar ='HSTREAM'
    CALL read_trip(nlisting,yfile_gw_frc_hd,yvar,zread)
    zhstream(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
    zhstream(:,:) = min(tp%XHC_BED(:,:),zhstream(:,:))
!
    yvar ='DRAIN'
    CALL read_trip(nlisting,yfile_gw_frc_hd,yvar,zread)
    zdrain(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end) !kg/m2/s
    zdrain(:,:) = max(0.0,zdrain(:,:)) * tpg%XAREA(:,:)      !kg/s
!
    tp%XWEFF   (:,:) = 0.0
! 
    DO jiter=1,2
       CALL gwf(tpg, &
                klon,klat,lwork,xtstep,xtstep,           &
                tpg%GMASK_GW,tp%XNUM_AQUI,zdrain,        &
                zlen,tp%XWIDTH,tp%XHC_BED,               &
                tp%XTOPO_RIV,tp%XTAUG,tpg%XAREA,         &
                tp%XTRANS,tp%XWEFF,                      &
                tp%XTABGW_F,tp%XTABGW_H,zhstream,        &
                tp%XHGROUND,zhg_old                      )
    ENDDO
!
    DEALLOCATE(zhstream)   
    DEALLOCATE(zdrain  )   
    DEALLOCATE(zhg_old )
!
  ENDIF         
!
  yvar ='HGROUND'   
  CALL write_trip(nlisting,yfile_prep,yvar,tpg%GMASK_GW,tp%XHGROUND,odouble=ldouble)
!
ENDIF
!
! * Calculate floodplains parameters
!
IF(lflood)THEN
!
  ALLOCATE(incell_bas(tpg%NBASMAX))
!
  incell_bas(:)=0
  DO jlat=1,ilat
     DO jlon=1,ilon
        jbas=tpg%NBASID(jlon,jlat)
        IF(jbas>0)THEN
          incell_bas(jbas)=incell_bas(jbas)+1
        ENDIF
     ENDDO
  ENDDO
!
  yvar ='NFLOOD'
  CALL read_trip(nlisting,yfile_read,yvar,zread)
  tp%XN_FLOOD(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
!
  ALLOCATE(zread3d(ilon_g,ilat_g,ndimtab))
  zread3d(:,:,:)= xundef
!  
  ALLOCATE(tp%XTAB_F (ilon,ilat,ndimtab))      
  ALLOCATE(tp%XTAB_H (ilon,ilat,ndimtab))      
  ALLOCATE(tp%XTAB_VF(ilon,ilat,ndimtab))  
!  
  tp%XTAB_F (:,:,:)= xundef
  tp%XTAB_H (:,:,:)= xundef
  tp%XTAB_VF(:,:,:)= xundef
!
  yvar ='TABF'
  CALL read_trip(nlisting,yfile_flood_read,yvar,zread3d(:,ilatf:ilat_g,:))
  tp%XTAB_F(:,:,:) = zread3d(ilon_deb:ilon_end,ilat_deb:ilat_end,:)
  yvar ='TABH'
  CALL read_trip(nlisting,yfile_flood_read,yvar,zread3d(:,ilatf:ilat_g,:))
  tp%XTAB_H(:,:,:) = zread3d(ilon_deb:ilon_end,ilat_deb:ilat_end,:)
  yvar ='TABVF'
  CALL read_trip(nlisting,yfile_flood_read,yvar,zread3d(:,ilatf:ilat_g,:))
  tp%XTAB_VF(:,:,:) = zread3d(ilon_deb:ilon_end,ilat_deb:ilat_end,:)  
!
  DEALLOCATE(zread3d)
!
  WHERE(tp%XHC_BED(:,:)==xundef)
        tp%XN_FLOOD(:,:)=xundef
  ENDWHERE
!
  WHERE(tp%XTAB_F(:,:,1)==xundef)
        tp%XN_FLOOD(:,:)=xundef
  ENDWHERE
!
  WHERE(tp%XTAB_H(:,:,1)==xundef)
        tp%XN_FLOOD(:,:)=xundef
  ENDWHERE
!
  lwork=(ilat*ilon>zmincell)
!  
  DO jlat=1,ilat
     DO jlon=1,ilon
        jbas=tpg%NBASID(jlon,jlat)
        IF(jbas==0)cycle
        
        IF(lwork.AND.incell_bas(jbas)<=zmincell)THEN
           tp%XN_FLOOD(jlon,jlat)=xundef
        ENDIF
     ENDDO
  ENDDO
!  
  tpg%GMASK_FLD(:,:)=tpg%GMASK_VEL(:,:)
  WHERE(tp%XN_FLOOD(:,:)==xundef)
        tpg%GMASK_FLD(:,:)=.false.
  ENDWHERE
!
  yvar ='NFLOOD'
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_FLD,tp%XN_FLOOD)
!
  yvar ='TABF'
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_FLD,tp%XTAB_F)
  yvar ='TABH'
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_FLD,tp%XTAB_H)  
  yvar ='TABVF'
  CALL write_trip(nlisting,yfile_param,yvar,tpg%GMASK_FLD,tp%XTAB_VF)
!
! Floodplains initial variables
!
  IF(lread_flood)THEN 
    yvar ='FLOOD_STO'
    CALL read_trip(nlisting,yfile_read_init,yvar,zread)
    tp%XFLOOD_STO(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
    yvar ='FFLOOD'
    CALL read_trip(nlisting,yfile_read_init,yvar,zread)
    tp%XFFLOOD(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
    yvar ='HFLOOD'
    CALL read_trip(nlisting,yfile_read_init,yvar,zread)
    tp%XHFLOOD(:,:) = zread(ilon_deb:ilon_end,ilat_deb:ilat_end)
  ELSE
    tp%XFLOOD_STO(:,:) = 0.0
    tp%XFFLOOD   (:,:) = 0.0
    tp%XHFLOOD   (:,:) = 0.0
  ENDIF
!
  yvar ='FLOOD_STO'
  CALL write_trip(nlisting,yfile_prep,yvar,tpg%GMASK,tp%XFLOOD_STO,odouble=ldouble)      
!
  yvar ='FFLOOD'
  CALL write_trip(nlisting,yfile_prep,yvar,tpg%GMASK,tp%XFFLOOD,odouble=ldouble)          
!
  yvar ='HFLOOD'
  CALL write_trip(nlisting,yfile_prep,yvar,tpg%GMASK,tp%XHFLOOD,odouble=ldouble)          
!
  DEALLOCATE(incell_bas)  
!
ENDIF
!
DEALLOCATE(zread)
!
!-------------------------------------------------------------------------------
! * Store current start date
!-------------------------------------------------------------------------------
!
 CALL write_trip_date(nlisting,yfile_prep,kyear,kmonth,kday,ptime)
!
!-------------------------------------------------------------------------------
! * END
!-------------------------------------------------------------------------------
!
WRITE(nlisting,*)''
WRITE(nlisting,*)'        END PREP TRIP            '
WRITE(nlisting,*)''
IF (lhook) CALL dr_hook('PREP_TRIP_RUN',1,zhook_handle)
!
!-------------------------------------------------------------------------------
END SUBROUTINE prep_trip_run