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 (USS, &
                              hprogram,hscheme,hsubroutine,hfilename,hfield)
!     #########################################################
!
!!**** *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_atm_sso_n, ONLY : surf_atm_sso_t
!
USE modd_pgd_grid,   ONLY : llatlonmask, xmeshlength
!
USE modd_arch, ONLY : little_endian_arch
!
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 yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
USE modi_refresh_pgdwork
USE modd_csts ,ONLY : xsurf_epsilon
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
!
TYPE(surf_atm_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.
!
!*    0.2    Declaration of local variables read in the data file head
!            ---------------------------------------------------------
!
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
INTEGER :: inbline                    ! number of latitude rows (number of lines
INTEGER :: inbcol                     ! number of longitude rows (number of columns)
REAL    :: znodata                    ! value below which data are not considered
!
!*    0.3    Declaration of local variables
!            ------------------------------
!
INTEGER :: iglb, iglbhdr              ! logical units
INTEGER :: iluout                     ! output listing logical unit
INTEGER :: ierr                       ! return codes
!
 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   ')
INTEGER           :: ibits            ! number of bits of a record in the
!                                     ! direct access file (16,32,64)
INTEGER :: jloop                      ! loop index
 CHARACTER(LEN=100):: ystring          ! string
 CHARACTER(LEN=88 ):: ystring1         ! part of string STRING
INTEGER           :: iindex           ! index of a character in string STRING1
!
REAL    :: zdlat                      ! latitude mesh in the data file
REAL    :: zdlon                      ! longitude mesh in the data file
INTEGER :: jline                      ! index of line
INTEGER :: jcol                       ! index of column
INTEGER, DIMENSION(2) :: icol1, icol2 ! limits of index of columns
INTEGER               :: iline1,iline2! limits of index of lines
INTEGER               :: icol         ! number of columns in mask domain
INTEGER               :: icolindex    ! column index in record
REAL, DIMENSION(:), POINTER :: zlat   ! latitude of data points
REAL, DIMENSION(:), POINTER :: zlon   ! longitude of data points
!
INTEGER :: jlon, jlat                 ! loop counters
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
REAL    :: zfact                      ! Factor integer to real
!
 CHARACTER,        DIMENSION(:), ALLOCATABLE :: ivalue8 ! value of a data point
 CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: ivalue16 ! value of a data point
 CHARACTER(LEN=4), DIMENSION(:), ALLOCATABLE :: ivalue32r ! value of a data point
INTEGER (KIND=4), DIMENSION(:), ALLOCATABLE :: ivalue32 ! value of a data point
INTEGER (KIND=8), DIMENSION(:), ALLOCATABLE :: ivalue64 ! value of a data point
REAL    (KIND=4), DIMENSION(:), ALLOCATABLE :: zvalue32 ! value of a data point
 CHARACTER(LEN=8), DIMENSION(:), ALLOCATABLE :: zvalue64 ! value of a data point
!
REAL, DIMENSION(:), ALLOCATABLE   :: zvalue             ! 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                           :: iwork              ! index of these data
LOGICAL                           :: gswap              ! T: swap has been done
!
!
INTEGER          :: ireclength        ! record length
INTEGER          :: irec              ! record number
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,9
  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
!
IF(ytype=='INTEGER')THEN
  IF(hfield=='CTI'.OR.hfield=='sand fraction'.OR. &
     hfield=='clay fraction'.OR.hfield=='organic carbon')THEN
    zfact=100.0
  ELSEIF (hfield=='water depth') THEN
    zfact=10.0
  ELSE
    zfact=1.0
  ENDIF
ELSE
  zfact=1.0
ENDIF
!
!----------------------------------------------------------------------------
!
!*    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)  
IF (ierr/=0) CALL abor1_sfx('READ_DIRECT: PB IN FILE HEADER')
!
!*    2.2    Closing of header
!            -----------------
!
 CALL close_namelist(hprogram,iglbhdr)
!
!*    2.3    Dimension of work arrays
!            ------------------------
!
ALLOCATE(zlat_work(inbcol))
ALLOCATE(zlon_work(inbcol))
ALLOCATE(zvalue_work(inbcol))
!
!----------------------------------------------------------------------------
!
!*    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))
ALLOCATE (ivalue8(inbcol))
ALLOCATE (ivalue16(inbcol))
ALLOCATE (ivalue32(inbcol))
ALLOCATE (ivalue32r(inbcol))
ALLOCATE (ivalue64(inbcol))
ALLOCATE (zvalue32(inbcol))
ALLOCATE (zvalue64(inbcol))
!
!*    7.2    Openning of direct access file
!            ------------------------------
!
 CALL open_file(hprogram,iglb,yfilename,'UNFORMATTED',           &
                 haction='READ',haccess='DIRECT',krecl=ireclength )  
!
!----------------------------------------------------------------------------
!
!*    4.     loop on mask meshes (lat)
!            -------------------
!
gswap = .false.
!
jlat = 0
!
DO 
!
  jlat = jlat + 1
  IF (jlat==361) EXIT
!
  IF ( .NOT. any(llatlonmask(:,jlat)) ) cycle
!
  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)
  IF ( .NOT. any(zlat(:)<zlatmax .AND. zlat(:)>=zlatmin) ) cycle
!
!----------------------------------------------------------------------------
!
!*    8.     Loop on lines
!            -------------
!
  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) ivalue8(:)
      zvalue(:)=ivalue8(:)
      ! 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
      READ(iglb,rec=irec) ivalue16(:)
      zvalue(:)=ivalue16(:)
      IF (      any(abs(zvalue)>15000)   ) THEN
        IF (gswap) CALL abor1_sfx('READ_DIRECT: SWAP ALREADY DONE, CANNOT BE REDONE')
        little_endian_arch = .NOT. little_endian_arch
        gswap = .true.
        WRITE(iluout,*) '*******************************************************************'
        WRITE(iluout,*) 'Architecture of the machine needs to swap LITTLE_ENDIAN_ARCH to ', &
                           little_endian_arch  
        WRITE(iluout,*) '*******************************************************************'
        jlat = 0
        CALL refresh_pgdwork
        EXIT   ! rereads the file
      END IF

    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) ivalue32r(:)
      zvalue(:)=ivalue32r(:)
      IF (      any(abs(zvalue)>0. .AND. abs(zvalue)<1.e-50) &
           .OR. any(abs(zvalue)>1.e20)                       ) THEN
        IF (gswap) CALL abor1_sfx('READ_DIRECT: SWAP ALREADY DONE, CANNOT BE REDONE')
        little_endian_arch = .NOT. little_endian_arch
        gswap = .true.
        WRITE(iluout,*) '*******************************************************************'
        WRITE(iluout,*) 'Architecture of the machine needs to swap LITTLE_ENDIAN_ARCH to ', &
                         little_endian_arch
        WRITE(iluout,*) '*******************************************************************'
        jlat = 0
        CALL refresh_pgdwork
        EXIT
      END IF                
    ELSE IF (ytype=='REAL   ' .AND. ibits==64) THEN
      READ(iglb,rec=irec) zvalue64(:)
      zvalue(:)=zvalue64(:)
      IF (      any(abs(zvalue)>0. .AND. abs(zvalue)<1.e-50) &
             .OR. any(abs(zvalue)>1.e20)                       ) THEN  
        IF (gswap) CALL abor1_sfx('READ_DIRECT: SWAP ALREADY DONE, CANNOT BE REDONE')
        little_endian_arch = .NOT. little_endian_arch
        gswap = .true.
        WRITE(iluout,*) '*******************************************************************'
        WRITE(iluout,*) 'Architecture of the machine needs to swap LITTLE_ENDIAN_ARCH to ', &
                           little_endian_arch  
        WRITE(iluout,*) '*******************************************************************'
        jlat = 0
        CALL refresh_pgdwork
        EXIT
      END IF
    ELSE
      CALL abor1_sfx('READ_DIRECT1: DATA TYPE NOT SUPPORTED')
    END IF
!
    IF(hfield=='CTI')THEN
      WHERE(zvalue(:)<0.0)zvalue(:)=znodata
    ENDIF
!
    WHERE(zvalue(:)/=znodata)zvalue(:)=zvalue(:)/zfact
!
!----------------------------------------------------------------------------
!
!*    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
!
!-------------------------------------------------------------------------------
!
!*   12.     Call to the adequate subroutine (point by point treatment)
!            ----------------------------------------------------------
!
          IF (iwork>0) &
            CALL pt_by_pt_treatment(uss, &
                                    iluout, zlat_work(1:iwork),zlon_work(1:iwork), &
                                    zvalue_work(1:iwork),                          &
                                    hsubroutine                                    )  
!
!-------------------------------------------------------------------------------
      END DO
!-------------------------------------------------------------------------------
    END DO
!-------------------------------------------------------------------------------
  END DO
!-------------------------------------------------------------------------------
END DO
!
!-------------------------------------------------------------------------------
!

DEALLOCATE(zlat)
DEALLOCATE(zlon)
!
DEALLOCATE(zlat_work  )
DEALLOCATE(zlon_work  )
DEALLOCATE(zvalue_work)
!
DEALLOCATE (zvalue)
DEALLOCATE (ivalue8 )
DEALLOCATE (ivalue16)
DEALLOCATE (ivalue32)
DEALLOCATE (ivalue32r)
DEALLOCATE (ivalue64)
DEALLOCATE (zvalue32)
DEALLOCATE (zvalue64)
!
 CALL close_file(hprogram,iglb)
IF (lhook) CALL dr_hook('READ_DIRECT',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE read_direct