read_direct_gauss.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_gauss (USS, &
                                    hprogram,hscheme,hsubroutine,hfilename,hfield)
!     #########################################################
!
!!**** *READ_DIRECT_GAUSS1* 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
!----------------------------------------------------------------------------
!
!*    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
!
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
!            ------------------------------
!
 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,        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    :: 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                    ! 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
REAL    :: zfact                      ! Factor integer to real
!
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 :: 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          :: 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 :: inblines, isize
INTEGER :: iwork              ! index of these data
INTEGER :: jloop, jlon, jlat, jline, jcol                    ! loop index
!
INTEGER, DIMENSION(2) :: icol1, icol2 ! limits of index of columns
!
INTEGER (KIND=4), DIMENSION(:), ALLOCATABLE :: ivalue32 ! value of a data point
INTEGER (KIND=8), DIMENSION(:), ALLOCATABLE :: ivalue64 ! value of a data point
!
LOGICAL                           :: gswap              ! T: swap has been done
!
REAL(KIND=JPRB) :: zhook_handle
!----------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_1',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
!
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_1',1,zhook_handle)
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_2',0,zhook_handle)
!----------------------------------------------------------------------------
!
!*    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_GAUSS: PB IN FILE HEADER')
!
!*    2.2    Closing of header
!            -----------------
!
 CALL close_namelist(hprogram,iglbhdr)
!
!*    2.3    Dimension of work arrays
!            ------------------------
!
isize = inbline*(inbcol/((zglblatmax-zglblatmin)*2.))
ALLOCATE(zlat_work(isize))
ALLOCATE(zlon_work(isize))
ALLOCATE(zvalue_work(isize))
!
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_2',1,zhook_handle)
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_3',0,zhook_handle)
!----------------------------------------------------------------------------
!
!*    3.     Adapt subgrid mesh to input file resolution
!            -------------------------------------------
!
IF (hsubroutine=='A_OROG') CALL ini_ssowork(xmeshlength,zdlat,zdlon)
!
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_3',1,zhook_handle)
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_4',0,zhook_handle)
!----------------------------------------------------------------------------
!
!*    7.     Openning of direct access file
!            ------------------------------
!
!*    7.1    Record length
!            -------------
!
ireclength = ibits/8 * inbcol
ALLOCATE (zvalue(inbcol))
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
!            ------------------------------
!
 CALL open_file(hprogram,iglb,yfilename,'UNFORMATTED',           &
                 haction='READ',haccess='DIRECT',krecl=ireclength ) 
!
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_4',1,zhook_handle)
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_5',0,zhook_handle)
!----------------------------------------------------------------------------
!
!*    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
  !            -------------
  !
  iwork=0
  DO jline = iline1,iline2
    !
    inblines = iline2 - iline1 + 1
    !----------------------------------------------------------------------------
    !
    !*   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
      READ(iglb,rec=irec) yvalue16(:)
      zvalue(:)=yvalue16(:)
      !
      IF (      any(abs(zvalue)>15000)   ) THEN
        IF (gswap) CALL abor1_sfx('READ_DIRECT_GAUSS: 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) yvalue32r(:)
      zvalue(:)=yvalue32r(:)
      !
      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_GAUSS: 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) yvalue64(:)
      zvalue(:)=yvalue64(:)
      !
      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_GAUSS: 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_DIRECT_GAUSS1: 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.e-10)/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
        !            ---------------
        !
        DO jcol=1,icol
          !
          !*   11.1    Recovers point value
          !            --------------------
          !
          icolindex = jcol+icol1(jloop)-1
          !
          !*   11.2    Test with respect to the 'no data' value
          !            ----------------------------------------
          !
          !*   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 ! JCOL
        !-------------------------------------------------------------------------------
      END DO !JLOOP
      !-------------------------------------------------------------------------------
    END DO !JLON
    !-------------------------------------------------------------------------------
  END DO ! JLINE
  !-------------------------------------------------------------------------------
  !
  !*   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, inblines, znodata                 )
!
!-------------------------------------------------------------------------------
END DO !JLAT
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_5',1,zhook_handle)
IF (lhook) CALL dr_hook('READ_DIRECT_GAUSS_6',0,zhook_handle)

DEALLOCATE(zlat)
DEALLOCATE(zlon)
!
DEALLOCATE(zlat_work  )
DEALLOCATE(zlon_work  )
DEALLOCATE(zvalue_work)
!
DEALLOCATE (zvalue)
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_GAUSS_6',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE read_direct_gauss