ol_read_atm_conf_netcdf.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 ol_read_atm_conf_netcdf (DTCO, U, HGRID, HSURF_FILETYPE, ODELAYEDSTART_NC, KDATESTOP,  &
                                    PDURATION, PTSTEP_FORC, KNI, KYEAR, KMONTH, KDAY, PTIME, &
                                    PLAT, PLON, PZS, PZREF, PUREF,KTIMESTARTINDEX         )  
!
!==================================================================
!!****  *OL_READ_ATM_CONF* - Initialization routine
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      F. Habets   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      Modified by P. Le Moigne (04/2005): cleaning and checking
!!      Modified by P. Le Moigne (04/2006): init_io_surf for nature
!!                  with GTMSK to read dimensions.
!!      Modified by Matthieu Lafaysse 2012-11-12
!==================================================================
!
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_surf_atm_n, ONLY : surf_atm_t
!
USE modd_type_date_surf
!
USE modd_surfex_mpi, ONLY : nrank, npio, ncomm, nproc, xtime_comm_read, xtime_npio_read
!
USE modi_set_surfex_filein
USE modi_get_luout
USE modi_init_io_surf_n
USE modi_read_surf
USE modi_end_io_surf_n
USE modi_get_size_full_n
USE modi_abor1_sfx
!
USE mode_dates_netcdf, ONLY : netcdf2date
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE netcdf
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include "mpif.h"
#endif
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(surf_atm_t), INTENT(INOUT) :: U
!
 CHARACTER(LEN=*), INTENT(IN)  :: HGRID
 CHARACTER(LEN=6), INTENT(IN)  :: HSURF_FILETYPE
LOGICAL, INTENT(IN)            :: ODELAYEDSTART_NC ! Allow the simulation to start from a different time step than the first record of a netcdf file
INTEGER,DIMENSION(4),INTENT(IN) :: KDATESTOP !Allow the simulation to end at a different time step than the last record of a netcdf file 
INTEGER,          INTENT(OUT) :: KNI
INTEGER,          INTENT(OUT) :: KYEAR, KMONTH, KDAY
REAL,             INTENT(OUT) :: PDURATION,PTSTEP_FORC
REAL,             INTENT(OUT) :: PTIME
REAL, DIMENSION(:),  POINTER  :: PLAT, PLON
REAL, DIMENSION(:),  POINTER  :: PZS 
REAL, DIMENSION(:),  POINTER  :: PZREF, PUREF
INTEGER,          INTENT(OUT) :: KTIMESTARTINDEX ! index from which we start reading FORCING.nc
!
 CHARACTER(LEN=100)           :: YUNITS
REAL, DIMENSION(:), POINTER   :: ZTIMEFILE
REAL,DIMENSION(:),ALLOCATABLE :: ZLAT1D, ZLON1D
REAL                          :: ZTIME
REAL                          :: ZFIRSTTIMEFILE
INTEGER                       :: IYEAR, IMONTH, IDAY
INTEGER                       :: IRET, INB_FORC
INTEGER                       :: INI, IDIM_FULL
INTEGER                       :: ILUOUT
INTEGER                       :: JINDEX
INTEGER                       :: INFOMPI
type(date_time)              :: ttime
type(date_time), DIMENSION(1) :: tz_firstdatefile
type(date_time), DIMENSION(:), ALLOCATABLE :: tz_datefile
INTEGER :: JPOINT, INLAT1D, INLON1D
LOGICAL :: GLONLAT1D
DOUBLE PRECISION   :: XTIME0
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!==================================================================
!
!*      0.    IO initialization
!
IF (lhook) CALL dr_hook('OL_READ_ATM_CONF_NETCDF',0,zhook_handle)
!
IF (nrank==npio) THEN
  !
#ifdef SFX_MPI
  xtime0 = mpi_wtime()
#endif
  !
  CALL get_luout(hsurf_filetype,iluout)
  !
  !*      1.    Read parameters from netcdf forcing file
  !
  yunits = ""
  CALL read_surf_dim_ol(yunits, inb_forc, ini, zfirsttimefile, iret,&
                        glonlat1d, inlon1d, inlat1d, ztimefile)
  !
#ifdef SFX_MPI
  xtime_npio_read = xtime_npio_read + (mpi_wtime() - xtime0)
#endif
  !
ENDIF
!
IF (nproc>1) THEN
#ifdef SFX_MPI
  xtime0 = mpi_wtime()
  CALL mpi_bcast(inb_forc,kind(inb_forc)/4,mpi_integer,npio,ncomm,infompi)
  CALL mpi_bcast(glonlat1d,1,mpi_logical,npio,ncomm,infompi)
  CALL mpi_bcast(inlon1d,kind(inlon1d)/4,mpi_integer,npio,ncomm,infompi)
  CALL mpi_bcast(inlat1d,kind(inlat1d)/4,mpi_integer,npio,ncomm,infompi)    
  xtime_comm_read = xtime_comm_read + (mpi_wtime() - xtime0)
#endif
ENDIF
!
 CALL read_surf(&
                'OFFLIN','FRC_TIME_STP'  ,ptstep_forc   ,iret)
!
pduration = ( inb_forc - 1 ) * ptstep_forc
!
!*      2.    Read full grid dimension and date
!
 CALL set_surfex_filein(hsurf_filetype,'PREP')
 CALL init_io_surf_n(dtco, u, hsurf_filetype,'FULL  ','SURF  ','READ ')  
!
 CALL read_surf(hsurf_filetype,'DIM_FULL',idim_full,iret)
 CALL read_surf(hsurf_filetype,'DTCUR',ttime,iret)
!
kyear  = ttime%TDATE%YEAR
kmonth = ttime%TDATE%MONTH
kday   = ttime%TDATE%DAY
ptime  = ttime%TIME
!
 CALL end_io_surf_n(hsurf_filetype)
!
!*      5.    Geographical initialization
!
 CALL get_size_full_n('OFFLIN ',idim_full,u%NSIZE_FULL,kni) 
!
ALLOCATE(plon(kni))
ALLOCATE(plat(kni))
ALLOCATE(pzs(kni))
ALLOCATE(pzref(kni))
ALLOCATE(puref(kni))
!
IF (glonlat1d) THEN
  ALLOCATE(zlat1d(inlat1d))
  ALLOCATE(zlon1d(inlon1d))
  ! we have to read PLAT and PLON due to parallel IO.
  ! however lat / lon have just INLAT1D and INLON1D dimension lengths in the netcdf file
  CALL read_surf('OFFLIN','LAT',plat,iret)
  CALL read_surf('OFFLIN','LON',plon,iret)
  zlat1d=plat(1:inlat1d)
  zlon1d=plon(1:inlon1d)
ELSE
  CALL read_surf('OFFLIN','LAT',plat,iret)
  CALL read_surf('OFFLIN','LON',plon,iret)
ENDIF
!
 CALL read_surf('OFFLIN','ZS',pzs,iret)
 CALL read_surf('OFFLIN','ZREF',pzref,iret)
 CALL read_surf('OFFLIN','UREF',puref,iret)
!
! Modif Matthieu Lafaysse
! If LON LAT REGULAR grid and 1 dimension LON LAT variables
! Affect 1d lon lat to total tables
IF (glonlat1d) THEN
  !
  IF ( inlat1d*inlon1d==idim_full ) THEN
    DO jpoint = 1,idim_full
      plat(jpoint) = zlat1d((jpoint-1)/inlon1d+1)
      plon(jpoint) = zlon1d(mod(jpoint-1,inlon1d)+1)
    END DO
    DEALLOCATE(zlat1d)
    DEALLOCATE(zlon1d)
  ELSE
    WRITE(iluout,*)' NUMBER OF GRID POINTS INCONSISTENCY: ',idim_full,'/',inlat1d*inlon1d
    CALL abor1_sfx('OL_READ_ATM_CONF_NETCDF: NUMBER OF GRID POINTS INCONSISTENCY')
  END IF
END IF
!
!*      6.    Check the consistency
!
IF (nrank == npio) THEN
  !
  IF (idim_full /= ini) THEN
    WRITE(iluout,*)' NUMBER OF GRID POINTS INCONSISTENCY: ',idim_full,'/',ini
    CALL abor1_sfx('OL_READ_ATM_CONF_NETCDF: NUMBER OF GRID POINTS INCONSISTENCY')
  ENDIF
  !
  IF ( odelayedstart_nc .OR. kdatestop(1)/=0 ) THEN
    !
    ! Convert time variable in dates
    ALLOCATE(tz_datefile(inb_forc))
    CALL netcdf2date(ztimefile, yunits, tz_datefile)
    !
  END IF
  DEALLOCATE(ztimefile)
  !
  IF ( odelayedstart_nc ) THEN
    !
    ktimestartindex = -1
    !
    ! Look for a date equal to the date prescribed in OPTIONS.NAM
    DO jindex = 1,inb_forc
      iyear  = tz_datefile(jindex)%TDATE%YEAR
      imonth = tz_datefile(jindex)%TDATE%MONTH
      iday   = tz_datefile(jindex)%TDATE%DAY
      ztime  = tz_datefile(jindex)%TIME* 3600. 
      
      IF ( kyear==iyear .AND. kmonth==imonth .AND. kday==iday .AND. ptime==ztime ) THEN
        ktimestartindex = jindex
        inb_forc        = inb_forc-ktimestartindex+1
        pduration       = ( inb_forc - 1 ) * ptstep_forc
        EXIT
      END IF
      
    END DO
    
    ! Error if the initial date is not found
    IF ( ktimestartindex==-1 ) THEN
      WRITE(iluout,*)'IN THE FORCING FILE, WE CAN NOT FIND THIS INITIAL DATE :',kyear,kmonth,kday,ptime
      CALL abor1_sfx('OL_READ_ATM_CONF_NETCDF: DATE INCONSISTENCY')
    END IF
    
  ELSE

    ktimestartindex = 1

    CALL netcdf2date((/zfirsttimefile/),yunits,tz_firstdatefile)
    iyear  = tz_firstdatefile(1)%TDATE%YEAR
    imonth = tz_firstdatefile(1)%TDATE%MONTH
    iday   = tz_firstdatefile(1)%TDATE%DAY
    ztime  = tz_firstdatefile(1)%TIME* 3600.
    !
    IF ( (kyear /= iyear) .OR. (kmonth /= imonth) .OR. (kday /= iday) ) THEN
      WRITE(iluout,*)' DATE INCONSISTENCY: ',kyear,kmonth,kday,'/',iyear,imonth,iday
      CALL abor1_sfx('OL_READ_ATM_CONF_NETCDF: DATE INCONSISTENCY')
    ENDIF
    !
    IF ( ptime /= ztime ) THEN
      WRITE(iluout,*)' TIME INCONSISTENCY: ',ptime,'/',ztime
      CALL abor1_sfx('OL_READ_ATM_CONF_NETCDF: TIME INCONSISTENCY')
    ENDIF
    !
  ENDIF
  !
  ! Look for the prescribed final date to modify PDURATION
  ! If we don't find it, ignore.
  IF ( kdatestop(1)/=0 ) THEN
    !
    DO jindex = ktimestartindex,inb_forc+ktimestartindex-1
      iyear  = tz_datefile(jindex)%TDATE%YEAR
      imonth = tz_datefile(jindex)%TDATE%MONTH
      iday   = tz_datefile(jindex)%TDATE%DAY
      ztime  = tz_datefile(jindex)%TIME* 3600.
      IF ( kdatestop(1)==iyear .AND. kdatestop(2)==imonth .AND. kdatestop(3)==iday &
                .AND. kdatestop(4)==int(ztime) ) THEN
        inb_forc  = jindex-ktimestartindex+1
        pduration = ( inb_forc - 1 ) * ptstep_forc
        EXIT
      END IF
    END DO
    !
  END IF
  !
  IF ( odelayedstart_nc .OR. kdatestop(1)/=0 ) THEN
    DEALLOCATE(tz_datefile)
  END IF  
  !
ENDIF
!
IF (nproc>1) THEN
#ifdef SFX_MPI       
  xtime0 = mpi_wtime()
  CALL mpi_bcast(ktimestartindex,kind(ktimestartindex)/4,mpi_integer,npio,ncomm,infompi)  
  CALL mpi_bcast(pduration,kind(pduration)/4,mpi_real,npio,ncomm,infompi)
  xtime_comm_read = xtime_comm_read + (mpi_wtime() - xtime0)
#endif
ENDIF
!
IF (lhook) CALL dr_hook('OL_READ_ATM_CONF_NETCDF',1,zhook_handle)
!==================================================================
CONTAINS
!
!     #############################################################
      SUBROUTINE read_surf_dim_ol(HUNITS,KSIZE,KNI,PFIRSTTIMEFILE,KRESP,&
                                  OLONLAT1D,KNLON1D,KNLAT1D,PTIMEFILE)
!     #############################################################
!
USE modi_ol_find_file_read
!
USE netcdf
!
IMPLICIT NONE
!
!*      0.1   Declarations of arguments
!
 CHARACTER(LEN=100), INTENT(OUT) :: HUNITS
INTEGER,            INTENT(OUT) :: KSIZE
INTEGER,            INTENT(OUT) :: KNI
REAL,               INTENT(OUT) :: PFIRSTTIMEFILE
INTEGER,            INTENT(OUT) :: KRESP    
LOGICAL,            INTENT(OUT) :: OLONLAT1D
INTEGER,            INTENT(OUT) :: KNLON1D,KNLAT1D
REAL,DIMENSION(:), POINTER, INTENT(OUT) :: PTIMEFILE
!
!*      0.2   Declarations of local variables
!
INTEGER :: IFILE_ID,IVAR_ID,INDIMS,JRET,JDIM,ITYPE
INTEGER,DIMENSION(2) :: IDIMIDS,IDIMLEN
 CHARACTER(20),DIMENSION(2) :: YDIMNAMES
INTEGER,DIMENSION(8) :: IRET

REAL*4, DIMENSION(:), ALLOCATABLE :: ZTIMEFILE4
REAL, DIMENSION(:), ALLOCATABLE   :: ZTIMEFILE
INTEGER, DIMENSION(:), ALLOCATABLE :: ITIMEFILE
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------

IF (lhook) CALL dr_hook('OL_READ_ATM_CONF_NETCDF:READ_SURF_DIM_OL',0,zhook_handle)
kresp=0

! 0. find filename
! -----------------
 CALL ol_find_file_read('time',ifile_id)
 
IF (ifile_id.NE.0) THEN
    
  ! 1. Find id of the variable
  !----------------------------
  iret(1)=nf90_inq_varid(ifile_id,'time',ivar_id)
  iret(2)=nf90_inquire_variable(ifile_id,ivar_id,ndims=indims)
  iret(3)=nf90_inquire_variable(ifile_id,ivar_id,dimids=idimids(1:1))
  iret(4)=nf90_inquire_dimension(ifile_id,idimids(1),len=ksize)
  iret(5)=nf90_get_att(ifile_id,ivar_id,'units',hunits)

  iret(6)=nf90_inquire_variable(ifile_id,ivar_id,xtype=itype)

  ALLOCATE(ptimefile(ksize))

  SELECT CASE (itype)
    CASE (nf90_double)
      ALLOCATE(ztimefile(ksize))
      iret(7)=nf90_get_var(ifile_id,ivar_id,ztimefile)
      pfirsttimefile=ztimefile(1)
      ptimefile(:)=ztimefile(:)
      DEALLOCATE(ztimefile)
    CASE (nf90_float)
      ALLOCATE(ztimefile4(ksize))
      iret(7)=REAL(nf90_get_var(ifile_id,ivar_id,ztimefile4))
      pfirsttimefile=ztimefile4(1)
      ptimefile(:)=ztimefile4(:)
      DEALLOCATE(ztimefile4)
    CASE (nf90_int)
      ALLOCATE(itimefile(ksize))
      iret(7)=REAL(nf90_get_var(ifile_id,ivar_id,itimefile))
      pfirsttimefile=itimefile(1)
      ptimefile(:)=itimefile(:)
      DEALLOCATE(itimefile)
    CASE DEFAULT
      CALL abor1_sfx('OL_READ_ATM_CONF_NETCDF: TYPE OF TIME VARIABLE NOT KNOWN')
  END SELECT
  
  ! 3. Check for errors
  !--------------------
  DO jret=1,7
    IF (ifile_id==0.OR.iret(jret).NE.nf90_noerr) kresp=1
  ENDDO

ENDIF

 CALL ol_find_file_read('LON',ifile_id)
!
olonlat1d = .false.
knlat1d = 0
knlon1d = 0
!
IF (ifile_id.NE.0) THEN

  iret(1)=nf90_inq_varid(ifile_id,'LON',ivar_id)
  iret(2)=nf90_inquire_variable(ifile_id,ivar_id,ndims=indims)
  iret(3)=nf90_inquire_variable(ifile_id,ivar_id,dimids=idimids(:))
  idimlen(:)=1.
  DO jdim=1,indims
    iret(4)=nf90_inquire_dimension(ifile_id,idimids(jdim),name=ydimnames(jdim),len=idimlen(jdim))
  ENDDO
!
  IF ( hgrid.EQ.'LONLAT REG' .AND. indims==1 ) THEN
    IF (trim(ydimnames(1))/="Number_of_points") THEN
      ! Modif Matthieu Lafaysse
      ! If LON LAT REGULAR grid and 1 dimension LON LAT variables
      ! total dimension is nlon*nlat
    
      olonlat1d = .true.
      CALL ol_find_file_read('LAT',ifile_id)
      iret(5) = nf90_inq_varid(ifile_id,'LAT',ivar_id)
      iret(6) = nf90_inquire_variable(ifile_id,ivar_id,ndims=indims)
      IF ( indims==1 ) THEN
        iret(7) = nf90_inquire_variable(ifile_id,ivar_id,dimids=idimids(2:2))
        iret(8) = nf90_inquire_dimension(ifile_id,idimids(2),len=idimlen(2))
        knlon1d = idimlen(1)
        knlat1d = idimlen(2)
      ELSE
        kresp = 1
      END IF
    END IF
  END IF
  kni = idimlen(1) * idimlen(2)

  DO jret=1,4
    IF (ifile_id==0.OR.iret(jret).NE.nf90_noerr) kresp=1
  ENDDO

ENDIF
!
IF (lhook) CALL dr_hook('OL_READ_ATM_CONF_NETCDF:READ_SURF_DIM_OL',1,zhook_handle)
!-------------------------------------------------------------------------------
END SUBROUTINE read_surf_dim_ol
!
END SUBROUTINE ol_read_atm_conf_netcdf