read_direct.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 read_direct (UG, U, USS, &
                              HPROGRAM,HSCHEME,HSUBROUTINE,HFILENAME,HFIELD,OMULTITYPE)
!     #########################################################
!
!!**** *READ_DIRECT1* reads a latlon file and call treatment subroutine
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!   
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson        Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    11/09/95
!!
!! V. Masson, March 2010     Optimization of some lat/lon boundaries computations
!!      J.Escobar     06/2013  for REAL4/8 add EPSILON management
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_surf_par, ONLY : xundef
USE modd_surfex_mpi, ONLY : nrank, nproc, npio
!
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_pgd_grid,   ONLY : llatlonmask, xmeshlength
!
USE modd_arch, ONLY : little_endian_arch
!
USE modd_data_cover_par, ONLY : jpcover, ntype
!
USE modi_get_luout
USE modi_open_namelist
USE modi_close_namelist
USE modi_open_file
USE modi_close_file
USE modi_readhead
USE modi_ini_ssowork
USE modi_pt_by_pt_treatment
USE mode_char2real
!
USE modi_uncompress_field
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
USE modd_pgdwork, ONLY : nsize_all, xall, nvalnbr, nvalcount, xvallist, &
                         catype, jpvalmax
USE modi_refresh_pgdwork
!
USE modd_csts ,ONLY : xsurf_epsilon
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(sso_t), INTENT(INOUT) :: USS
!
 CHARACTER(LEN=6),  INTENT(IN) :: HPROGRAM      ! Type of program
 CHARACTER(LEN=6),  INTENT(IN) :: HSCHEME       ! Scheme treated
 CHARACTER(LEN=6),  INTENT(IN) :: HSUBROUTINE   ! Name of the subroutine to call
 CHARACTER(LEN=28), INTENT(IN) :: HFILENAME     ! Name of the field file.
 CHARACTER(LEN=20), INTENT(IN) :: HFIELD        ! Name of the field.
LOGICAL, OPTIONAL, INTENT(IN) :: OMULTITYPE
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
 CHARACTER(LEN=28) :: YFILENAME        ! Name of the field file without header
 CHARACTER(LEN=28) :: YFILEHDR         ! Name of the field file header
!
 CHARACTER(LEN=7)  :: YTYPE            ! type of numerical field stored in the
!                                     ! direct access file ('INTEGER','REAL   ')
 CHARACTER(LEN=100):: YSTRING          ! string
 CHARACTER(LEN=88 ):: YSTRING1         ! part of string STRING
!
 CHARACTER(LEN=2), DIMENSION(1) :: YCPT16 ! value of a data point 
 CHARACTER(LEN=4), DIMENSION(:), ALLOCATABLE :: YCPT32
!
 CHARACTER,        DIMENSION(:), ALLOCATABLE :: YVALUE8 ! value of a data point
 CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: YVALUE16 ! value of a data point
 CHARACTER(LEN=4), DIMENSION(:), ALLOCATABLE :: YVALUE32R ! value of a data point
 CHARACTER(LEN=8), DIMENSION(:), ALLOCATABLE :: YVALUE64 ! value of a data point
!
REAL, DIMENSION(1) :: ZCPT
!
REAL    :: ZGLBLATMIN                 ! minimum latitude of data box in the file
REAL    :: ZGLBLONMIN                 ! minimum longitude of data box in the file
REAL    :: ZGLBLATMAX                 ! maximum latitude of data box in the file
REAL    :: ZGLBLONMAX                 ! maximum longitude of data box in the file
REAL    :: ZNODATA, ZNODATA2          ! value below which data are not considered
REAL    :: ZDLAT                      ! latitude mesh in the data file
REAL    :: ZDLON                      ! longitude mesh in the data file
REAL    :: ZLONMIN                    ! minimum longitude of mask mesh
REAL    :: ZLONMAX                    ! maximum longitude of mask mesh
REAL    :: ZLATMIN                    ! minimum latitude of mask mesh
REAL    :: ZLATMAX                    ! maximum latitude of mask mesh
REAL    :: ZSHIFT                     ! shift on longitudes
INTEGER :: IFACT              ! Factor integer to real
INTEGER(KIND=2) :: INODATA, INODATA2
!
REAL, DIMENSION(:), ALLOCATABLE :: ZVALUE
REAL, DIMENSION(:), POINTER :: ZLAT   ! latitude of data points
REAL, DIMENSION(:), POINTER :: ZLON   ! longitude of data points
REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: ZVALUE32 ! value of a data point
REAL,         DIMENSION(:), ALLOCATABLE :: ZINTER     ! value of a record of data points
REAL,         DIMENSION(:), ALLOCATABLE :: ZVALUE_WORK        ! value of a valid data points 
REAL,         DIMENSION(:), ALLOCATABLE :: ZLAT_WORK          ! latitude  of a valid data points 
REAL,         DIMENSION(:), ALLOCATABLE :: ZLON_WORK          ! longitude of a valid data points 
!
INTEGER          :: IINDEX           ! index of a character in string STRING1
INTEGER          :: IBITS             ! number of bits of a record in the
                                      ! direct access file (16,32,64)
INTEGER          :: IRECLENGTH        ! record length
INTEGER          :: IREC              ! record number
INTEGER          :: IGLB, IGLBHDR     ! logical units
INTEGER          :: ILUOUT            ! output listing logical unit
INTEGER          :: IERR              ! return codes
!
INTEGER :: INBLINE                    ! number of latitude rows (number of lines
INTEGER :: INBCOL                     ! number of longitude rows (number of columns)
INTEGER :: ILINE1,ILINE2              ! limits of index of lines
INTEGER :: ICOL                       ! number of columns in mask domain
INTEGER :: ICOLINDEX                  ! column index in record
INTEGER :: ISIZE
INTEGER :: IWORK, IDEB, IPAS         ! index of these data
INTEGER :: JLOOP, JLON, JLAT, JLINE, JCOL, JL, ICPT, INB, JTYPE
INTEGER :: ILINE_COMPRESS
INTEGER :: INB_LINE_READ
INTEGER(KIND=8) :: IPOS
!
INTEGER, DIMENSION(360) :: IMASK
INTEGER, DIMENSION(2) :: ICOL1, ICOL2 ! limits of index of columns
!
INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ICPT0          ! loop index
INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: IVALUE32 ! value of a data point
INTEGER(KIND=8), DIMENSION(:), ALLOCATABLE :: IVALUE64 ! value of a data point
!
 CHARACTER(LEN=6) :: YACCESS
LOGICAL           :: GSWAP              ! T: swap has been done
LOGICAL :: GMULTITYPE, GCOMPRESS
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!----------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('READ_DIRECT',0,zhook_handle)
 CALL get_luout(hprogram,iluout)
!
!*    1.     Openning of global field
!            ------------------------
!
!*    1.1    Logical unit attributions
!            -------------------------
!
yfilename=adjustl(adjustr(hfilename)//'.dir')
yfilehdr =adjustl(adjustr(hfilename)//'.hdr')
!
!*    1.2    Openning of header
!            ------------------
!
 CALL open_namelist(hprogram,iglbhdr,yfilehdr)
!
!*    1.3    Reading in header of direct access characteristics
!            --------------------------------------------------
!
DO jloop=1,11
  READ(iglbhdr,'(A100)') ystring
  IF (ystring(1:10)=='recordtype') EXIT
END DO
!
rewind(iglbhdr)
!
ystring1=ystring(12:100)
!
!* string analysis
!
iindex=index(ystring1,'n')  ! n for integer
IF (iindex/=0) THEN
  ytype='INTEGER'
ELSE
  ytype='REAL   '
END IF
iindex=index(ystring1,'8')
IF (iindex/=0) ibits=8
iindex=index(ystring1,'1')
IF (iindex/=0) ibits=16
iindex=index(ystring1,'3')
IF (iindex/=0) ibits=32
iindex=index(ystring1,'4')
IF (iindex/=0) ibits=64
!
!----------------------------------------------------------------------------
!
!*    2.     Reading of the global field
!            ---------------------------
!
!*    2.1    Head of data file
!            -----------------
!
 CALL readhead(iglbhdr,zglblatmin,zglblatmax,zglblonmin,zglblonmax, &
               inbline,inbcol,znodata,zdlat,zdlon,zlat,zlon,ierr,ifact,&
               gcompress)  
IF (ierr/=0) CALL abor1_sfx('READ_DIRECT: PB IN FILE HEADER')
!
IF (gcompress .AND. (ytype/='INTEGER' .OR. ibits/=16)) &
  CALL abor1_sfx('READ_DIRECT: COMPRESSED FILES ARE POSSIBLE ONLY WITH INTEGER 16 BYTES FOR THE MOMENT')
!
gmultitype = .false.
IF (PRESENT(omultitype)) gmultitype = omultitype
!
IF (gmultitype) THEN
  DEALLOCATE(nsize_all)
  ALLOCATE(nsize_all(u%NDIM_FULL,sum(ntype)))  
  nsize_all(:,:) = 0
  IF (catype=='MAJ') THEN
    DEALLOCATE(nvalnbr,nvalcount,xvallist)
    ALLOCATE(nvalnbr(u%NDIM_FULL,sum(ntype)))
    ALLOCATE(nvalcount(u%NDIM_FULL,jpvalmax,sum(ntype)))
    ALLOCATE(xvallist(u%NDIM_FULL,jpvalmax,sum(ntype)))    
    nvalnbr(:,:)   = 0
    nvalcount(:,:,:) = 0
    xvallist(:,:,:) = xundef    
  ELSE
    DEALLOCATE(xall)
    ALLOCATE(xall(u%NDIM_FULL,sum(ntype),1))
    xall(:,:,:) = 0.
  ENDIF
ENDIF
!
IF(ytype=='INTEGER')THEN
  IF(hfield(1:3)=='CTI'.OR.hfield=='sand fraction'.OR.hfield=='clay fraction'.OR.&
     hfield=='organic carbon'.OR.hfield(1:4)=='SAND'.OR. hfield(1:4)=='CLAY'.OR.hfield(1:3)=='SOC')THEN
    ifact=100
  ELSEIF (hfield=='water depth') THEN
    ifact=10
  ENDIF
ENDIF
!
!*    2.2    Closing of header
!            -----------------
!
 CALL close_namelist(hprogram,iglbhdr)
!
!*    2.3    Dimension of work arrays
!            ------------------------
!
! ires c'est le nombre de lignes qu'on lit dans un demi degré (multiple de 60)
inb_line_read = inbline / ((zglblatmax-zglblatmin)*2.)
IF (inb_line_read>60) inb_line_read = max(inb_line_read/3,60)
! on lit toujours 60 lignes d'un coup 
isize = inb_line_read * inbcol
!
ALLOCATE(zlat_work(isize))
ALLOCATE(zlon_work(isize))
ALLOCATE(zvalue_work(isize))
IF (gcompress.OR.gmultitype) THEN
  ALLOCATE (zinter(isize))
  zinter(:) = 0.
ENDIF
!
!----------------------------------------------------------------------------
!
!*    3.     Adapt subgrid mesh to input file resolution
!            -------------------------------------------
!
IF (hsubroutine=='A_OROG') CALL ini_ssowork(xmeshlength,zdlat,zdlon)
!
!----------------------------------------------------------------------------
!
!*    7.     Openning of direct access file
!            ------------------------------
!
!*    7.1    Record length
!            -------------
!
ireclength = ibits/8 * inbcol
ALLOCATE (zvalue(inbcol))
zvalue(:) = 0.
!
IF (gcompress) THEN
  ALLOCATE(ycpt32(inbline))
  ALLOCATE(icpt0(inbline))
ENDIF
!
IF (ytype=='INTEGER' .AND. ibits== 8) THEN  
  ALLOCATE (yvalue8(inbcol))
ELSEIF (ytype=='INTEGER' .AND. ibits==16) THEN
  ALLOCATE (yvalue16(inbcol))
ELSEIF (ytype=='INTEGER' .AND. ibits==32) THEN
  ALLOCATE (ivalue32(inbcol))
  ELSEIF (ytype=='INTEGER' .AND. ibits==64) THEN
  ALLOCATE (ivalue64(inbcol))
ELSEIF (ytype=='REAL   ' .AND. ibits==32) THEN
  ALLOCATE (yvalue32r(inbcol))
ELSEIF (ytype=='REAL   ' .AND. ibits==64) THEN
  ALLOCATE (yvalue64(inbcol))
ENDIF
!
!*    7.2    Openning of direct access file
!            ------------------------------
!
yaccess = 'DIRECT'
IF (gcompress) THEN
  yaccess='STREAM'
  little_endian_arch = .false.
ENDIF
!
 CALL open_file(hprogram,iglb,yfilename,'UNFORMATTED',           &
                 haction='READ',haccess=yaccess,krecl=ireclength ) 
!
! we read numbers of elements by line of the grid at the beginning
IF (gcompress) THEN
  READ(iglb) ycpt32
  icpt0(:) = transfer(ycpt32(:),1_4,inbline)
ENDIF
!
!----------------------------------------------------------------------------
!
!*    4.     loop on mask meshes (lat)
!            -------------------
!
imask(:) = 0
icpt = 0
DO jlat = 1,360
  IF ( .NOT. any(llatlonmask(:,jlat)) ) cycle
  zlatmin = (jlat-180)/2. - 0.5
  zlatmax = (jlat-180)/2.
  IF ( .NOT. any(zlat(:)<zlatmax .AND. zlat(:)>=zlatmin) ) cycle
  icpt = icpt + 1
  imask(icpt) = jlat
ENDDO
!
!INB: number of lat to be read
inb = icpt
!
!IPAS: number of lat to be read for each task
ipas = ceiling(inb*1./nproc)
!
gswap = .false.
!
!first lat read for this task
ideb = ipas*nrank
!
icpt = 0
!
jl = ipas + 1
!
IF (gcompress) iline_compress = 1
!
inodata = znodata
inodata2 = ishftc(inodata,8)
znodata2 = inodata2
!
DO 
  !
  !the file is read from the top to the bottom (quicker)
  jl = jl - 1
  IF (jl==0) EXIT
  !
  IF (ideb+jl>inb) cycle
  !
  !lat read by this task for this loop index JL
  jlat = imask(ideb+jl)
  !
  zlatmin = (jlat-180)/2. - 0.5
  zlatmax = (jlat-180)/2.
  !
  !----------------------------------------------------------------------------
  !
  !*    5.     index limits on latitude
  !            ------------------------
  !
  iline1=max(min(int((zglblatmax-zdlat/2.-zlatmax)/zdlat+1.),inbline),0)+1
  iline2=max(min(int((zglblatmax-zdlat/2.-zlatmin)/zdlat+1.),inbline),0)
  !
  !----------------------------------------------------------------------------
  !
  !*    8.     Loop on lines
  !            -------------
  !
  ! first IPOS for this task is the first information plus the
  ! number of elements by lines before the first to read
  IF (gcompress.AND.(jl==ipas.OR.iline_compress<iline1)) THEN
    ipos = 0
    IF (iline1>1) THEN
      DO jloop=1,iline1-1
        ipos = ipos + icpt0(jloop)
      ENDDO
    ENDIF
    ipos = ipos*2 + 1 + inbline*4
    iline_compress = iline1
  ELSE
    ipos = 0
  ENDIF
  !
  DO jline = iline1,iline2
    !
    !----------------------------------------------------------------------------
    !
    !*   10.     Reading in the direct access file
    !            ---------------------------------
    !
    !*   10.1    Record number
    !            -------------
    !
    irec = jline
    ! 
    !*   10.2    Reading the correct data type and conversion into real
    !            ------------------------------------------------------
    !
    IF      (ytype=='INTEGER' .AND. ibits== 8) THEN
      READ(iglb,rec=irec) yvalue8(:)
      zvalue(:)=yvalue8(:)
      ! negative values are shifted to positive values according to binary coding
      WHERE (zvalue(:)<0.) zvalue(:) = nint(256.+zvalue(:))
      !
    ELSE IF (ytype=='INTEGER' .AND. ibits==16) THEN
      !
      IF (gcompress) THEN
        IF (ipos/=0) THEN
          READ(iglb,pos=ipos) yvalue16(1:icpt0(jline))
          ipos = 0
        ELSE
          READ(iglb) yvalue16(1:icpt0(jline))
        ENDIF
        zvalue(1:icpt0(jline))=yvalue16(1:icpt0(jline))
        iline_compress = iline_compress + 1
      ELSE
        READ(iglb,rec=irec) yvalue16(:)
        zvalue(:)=yvalue16(:)
      ENDIF
      !
      IF (icpt==0.AND..NOT.gcompress) THEN 
        IF ( (hfield(1:5)=="COVER" .AND. (any(zvalue>jpcover.AND.zvalue/=znodata) .OR. &
                        any(zvalue<0..AND.zvalue/=znodata) .OR. all(zvalue==256.)) ) .OR. & 
             (znodata/=0 .AND. (all(zvalue==znodata2))) .OR. &
            ((hfield(1:4)=="SAND" .OR. hfield(1:4)=="CLAY") .AND. &
                (any(zvalue>100..AND.zvalue/=znodata) .OR. any(zvalue<0..AND.zvalue/=znodata)) ) .OR. &
             (hfield(1:3)=="SOC" .AND. (any(zvalue>15000..AND.zvalue/=znodata) .OR. any(zvalue<0..AND.zvalue/=znodata)) )  .OR. &
             ((hfield(1:5)/="COVER" .AND. hfield(1:4)/="SAND" .AND. hfield(1:4)/="CLAY" .AND. &
              hfield(1:3)/="SOC" .AND. any(zvalue>15000..AND.zvalue/=znodata) ) ) ) THEN
          icpt = icpt + 1
          IF (gswap) CALL abor1_sfx('READ_DIRECT: SWAP ALREADY DONE, CANNOT BE REDONE')
          little_endian_arch = .NOT. little_endian_arch
          gswap = .true.
          IF (nrank==npio) THEN
            WRITE(iluout,*) '*******************************************************************'
            WRITE(iluout,*) 'Architecture of the machine needs to swap LITTLE_ENDIAN_ARCH to ', &
                             little_endian_arch  
            WRITE(iluout,*) '*******************************************************************'
          ENDIF
          jl = ipas + 1 !back to first lat
          IF (hfield(1:5)=="COVER") u%LCOVER(:) = .false.
          CALL refresh_pgdwork(hsubroutine)
          EXIT   ! rereads the file
        ENDIF
      ENDIF
      !
      !
    ELSE IF (ytype=='INTEGER' .AND. ibits==32) THEN
      READ(iglb,rec=irec) ivalue32(:)
      zvalue(:)=ivalue32(:)
      !
    ELSE IF (ytype=='INTEGER' .AND. ibits==64) THEN
      READ(iglb,rec=irec) ivalue64(:)
      zvalue(:)=ivalue64(:)
      !
    ELSE IF (ytype=='REAL   ' .AND. ibits==32) THEN
      READ(iglb,rec=irec) yvalue32r(:)
      zvalue(:)=yvalue32r(:)
      !
      IF (icpt==0) THEN      
        IF (      any(abs(zvalue)>0. .AND. abs(zvalue)<1.e-50) &
             .OR. any(abs(zvalue)>1.e20)                       ) THEN
          icpt = icpt + 1
          IF (gswap) CALL abor1_sfx('READ_DIRECT: SWAP ALREADY DONE, CANNOT BE REDONE')
          little_endian_arch = .NOT. little_endian_arch
          gswap = .true.
          IF (nrank==npio) THEN
            WRITE(iluout,*) '*******************************************************************'
            WRITE(iluout,*) 'Architecture of the machine needs to swap LITTLE_ENDIAN_ARCH to ', &
                             little_endian_arch
            WRITE(iluout,*) '*******************************************************************'
          ENDIF
          jl = ipas + 1
          CALL refresh_pgdwork(hsubroutine)
          EXIT
        ENDIF
      END IF      
      !
    ELSE IF (ytype=='REAL   ' .AND. ibits==64) THEN
      READ(iglb,rec=irec) yvalue64(:)
      zvalue(:)=yvalue64(:)
      !
      IF (icpt==0) THEN      
        IF (      any(abs(zvalue)>0. .AND. abs(zvalue)<1.e-50) &
               .OR. any(abs(zvalue)>1.e20)                       ) THEN  
          icpt = icpt + 1
          IF (gswap) CALL abor1_sfx('READ_DIRECT: SWAP ALREADY DONE, CANNOT BE REDONE')
          little_endian_arch = .NOT. little_endian_arch
          gswap = .true.
          IF (nrank==npio) THEN
            WRITE(iluout,*) '*******************************************************************'
            WRITE(iluout,*) 'Architecture of the machine needs to swap LITTLE_ENDIAN_ARCH to ', &
                             little_endian_arch  
            WRITE(iluout,*) '*******************************************************************'
          ENDIF
          jl = ipas + 1
          CALL refresh_pgdwork(hsubroutine)
          EXIT
        ENDIF
      ENDIF
    ELSE
      CALL abor1_sfx('READ_DIRECT1: DATA TYPE NOT SUPPORTED')
    END IF
    !
    IF(hfield=='CTI')THEN
      WHERE(zvalue(:)<0.0) zvalue(:)=znodata
    ENDIF
    !
    IF (gcompress) THEN
      WHERE (zvalue(1:icpt0(jline))<0.) zvalue(1:icpt0(jline)) = nint(32768*2.+zvalue(1:icpt0(jline)))
      zinter(1:icpt0(jline)) = zvalue(1:icpt0(jline))
      CALL uncompress_field(inbcol,4000.,zinter(1:icpt0(jline)),zvalue(:))
    ENDIF
    !
    !----------------------------------------------------------------------------
    !
    !*    4.     loop on mask meshes (lon)
    !            -------------------
    !
    DO jlon=1,720
      !
      IF (.NOT. llatlonmask(jlon,jlat)) cycle
      !
      zlonmin =  jlon     /2. - 0.5
      zlonmax =  jlon     /2.
      !
      !----------------------------------------------------------------------------
      !
      !*    5.     limits on longitude
      !            -------------------
      !
      !*    5.1    left domain border is set just higher than the global field min. longitude
      !            -----------------------------------------------------------------
      !
      zshift = 360. * nint((zlonmin-zglblonmin-180.*(1-xsurf_epsilon))/360.)
      !
      zglblonmin = zglblonmin + zshift
      zglblonmax = zglblonmax + zshift
      !
      !*    5.2    index limits on longitude
      !            -------------------------
      !
      icol1(1)=max(min(int((zlonmin-zglblonmin-zdlon/2.)/zdlon+1.),inbcol),0)+1
      icol2(1)=max(min(int((zlonmax-zglblonmin-zdlon/2.)/zdlon+1.),inbcol),0)
      !
      !* Does right domain border goes outside the global field longitude range?
      !* Does it then goes into the global field domain by the other side?
      !* Then a second part of the global field domain must be considered
      !
      icol1(2)=1
      icol2(2)=max(min(int((zlonmax-zglblonmin-zdlon/2.-360.)/zdlon+1.),inbcol),0)
      !
      !----------------------------------------------------------------------------
      !
      !*    6.     Loop on longitude limits
      !            ------------------------
      !
      DO jloop=1,2
        !
        icol = icol2(jloop) - icol1(jloop) + 1
        !
        IF (icol<1) cycle
        !
        !----------------------------------------------------------------------------
        !
        !*   11.     Loop on columns
        !            ---------------
        !
        iwork=0
        !
        DO jcol=1,icol
          !
          !*   11.1    Recovers point value
          !            --------------------
          !
          icolindex = jcol+icol1(jloop)-1
          !
          !*   11.2    Test with respect to the 'no data' value
          !            ----------------------------------------
          !
          IF (abs(zvalue(icolindex)-znodata)<=1.e-10) cycle
          !
          !*   11.3    copy of the correct values in a work array
          !            ------------------------------------------
          !
          iwork = iwork + 1
          zlat_work(iwork) = zlat(jline)
          zlon_work(iwork) = zlon(icolindex)
          zvalue_work(iwork) = zvalue(icolindex)
          ! 
        END DO
        !-------------------------------------------------------------------------------
        !
        IF (.NOT.gmultitype.AND.ifact/=1) THEN
          WHERE(zvalue_work(1:iwork)/=znodata) 
            zvalue_work(1:iwork)=zvalue_work(1:iwork)/float(ifact)
          END WHERE
        ENDIF
        !
        !*   12.     Call to the adequate subroutine (point by point treatment)
        !            ----------------------------------------------------------
        !
        IF (iwork>0) &
          CALL pt_by_pt_treatment(ug, u, uss, &
                                  iluout, zlat_work(1:iwork),zlon_work(1:iwork), &
                                  zvalue_work(1:iwork),                          &
                                  hsubroutine, omultitype=gmultitype, kfact=ifact)  
!
!-------------------------------------------------------------------------------
      END DO
!-------------------------------------------------------------------------------
    END DO
!-------------------------------------------------------------------------------
  END DO
!-------------------------------------------------------------------------------
END DO
!
!-------------------------------------------------------------------------------
!

DEALLOCATE(zlat)
DEALLOCATE(zlon)
!
DEALLOCATE(zlat_work  )
DEALLOCATE(zlon_work  )
DEALLOCATE(zvalue_work)
!
DEALLOCATE (zvalue)
IF (ALLOCATED(zinter)) DEALLOCATE (zinter)
IF (ALLOCATED(yvalue8)) DEALLOCATE (yvalue8 )
IF (ALLOCATED(yvalue16)) DEALLOCATE (yvalue16)
IF (ALLOCATED(yvalue32r)) DEALLOCATE (yvalue32r)
IF (ALLOCATED(yvalue64)) DEALLOCATE (yvalue64)
IF (ALLOCATED(ivalue32)) DEALLOCATE (ivalue32)
IF (ALLOCATED(ivalue64)) DEALLOCATE (ivalue64)
IF (ALLOCATED(zvalue32)) DEALLOCATE (zvalue32)
!
 CALL close_file(hprogram,iglb)
IF (lhook) CALL dr_hook('READ_DIRECT',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE read_direct