MODULE ReadNETCDF
 IMPLICIT NONE
 
! This is part of the netCDF package.
! Copyright 2006 University Corporation for Atmospheric Research/Unidata.
! See COPYRIGHT file for conditions of use.

! This is an example which reads some 4D pressure and
! temperatures. The data file read by this program is produced by
! the companion program pres_temp_4D_wr.f90. It is intended to
! illustrate the use of the netCDF Fortran 90 API.

! This program is part of the netCDF tutorial:
! http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-tutorial

! Full documentation of the netCDF Fortran 90 API can be found at:
! http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-f90

! $Id: pres_temp_4D_rd.f90,v 1.9 2010/04/06 19:32:09 ed Exp $
PRIVATE

  ! We are reading 4D data, a 12 x 6 x 2 lon-lat-lvl grid, with 2
  ! timesteps of data.
  integer, parameter :: nMaxFile   =   6
  integer, parameter :: nMaxFile2  =   7

  integer, parameter :: nvar  =   4
  integer, parameter :: NDIMS =   3
  integer,public, parameter :: NTIME = 73700
  integer, parameter :: NLATS =  96
  integer, parameter :: NLONS = 192
  integer ::  lon_dimid
  integer ::  lat_dimid
  integer ::  tim_dimid

  !global attribute
  CHARACTER (len = *), parameter :: history     = "history" 
  CHARACTER (len = *), parameter :: institution = "CPTEC" 
  CHARACTER (len = *), parameter :: source      = "source" 
  CHARACTER (len = *), parameter :: Conventions = "Conventions" 

  INTEGER, PARAMETER            :: MAX_ATT_LEN      = 800
  CHARACTER (len = MAX_ATT_LEN) :: gbl_history      != "NASA LaRC/GEWEX SRB" ;
  CHARACTER (len = MAX_ATT_LEN) :: gbl_Conventions  != "CF-1.1" ; 

  ! In addition to the latitude and longitude dimensions, we will also
  ! create latitude and longitude variables which will hold the actual
  ! latitudes and longitudes. Since they hold data about the
  ! coordinate system, the netCDF term for these is: "coordinate
  ! variables."
  REAL(KIND=8) :: lats     (  NLATS)
  REAL(KIND=8) :: lons     (  NLONS)
  REAL(KIND=8) :: timer    (  NTIME)

  ! We recommend that each variable carry a "units" attribute.
  TYPE :: attribute
     character (len =  5) :: UNITS              ="units"
     character (len = 10) :: FillValue          ="_FillValue"
     character (len = 10) :: add_offset         ="add_offset"
     character (len = 12) :: scale_factor       ="scale_factor"
     character (len = 13) :: standard_name      ="standard_name" 
     character (len =  9) :: long_name          ="long_name" 
     character (len = 12) :: actual_range       ="actual_range"
     character (len =  7) :: delta_t            ="delta_t"
     character(LEN=400)   :: var_units          != "W m-2" ;
     INTEGER(KIND=2)      :: var_FillValue      != -100s ;
     REAL                 :: lat_actual_range(2)! = -89.5f, 89.5f ;
     REAL                 :: lon_actual_range(2)! = -89.5f, 89.5f ;
     REAL                 :: var_add_offset     != 0.f ;
     REAL                 :: var_scale_factor   != 0.1f ;
     character(LEN=400)   :: var_standard_name  != "toa_incoming_shortwave_flux" ;
     character(LEN=400)   :: var_long_name      != "All-Sky Shortwave TOA Downward Flux" ;
     character(LEN=400)   :: time_delta_t       != "CF-1.1" ; 
     REAL(KIND=8)         :: time_actual_range(NTIME)       ! = -89.5f, 89.5f ;
     REAL(KIND=8)         :: time_actual_range_out(3*NTIME) ! = -89.5f, 89.5f ;
  END TYPE attribute
  TYPE(attribute), ALLOCATABLE :: var(:,:)
  ! The start and count arrays will tell the netCDF library where to
  ! read our data.
  integer :: start     (NDIMS)
  integer :: start2    (NDIMS)
  integer :: start_out (NDIMS)
  INTEGER :: count     (NDIMS)
  integer :: dimids    (NDIMS)
  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: var_c
  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: var_p
  ! We will read surface temperature and pressure fields. In netCDF
  ! terminology these are called "variables."
  CHARACTER (len = 20),PUBLIC, PARAMETER :: VAR_NAME(1:nvar,1:nMaxFile)=RESHAPE((/&
      "lon          ","lat          ","time         ","precip       ",&
      "lon          ","lat          ","time         ","temp         ",&
      "lon          ","lat          ","time         ","rhum         ",&
      "lon          ","lat          ","time         ","windspeed    ",&
      "lon          ","lat          ","time         ","swrad        ",&
      "lon          ","lat          ","time         ","lwrad        "/),(/nvar,nMaxFile/))
 
  CHARACTER (len = 20),PUBLIC, PARAMETER :: VAR_NAME2(1:nvar,1:nMaxFile2)=RESHAPE((/&
      "lon          ","lat          ","time         ","pres         ",&
      "lon          ","lat          ","time         ","tems         ",&
      "lon          ","lat          ","time         ","umes         ",&
      "lon          ","lat          ","time         ","wind         ",&
      "lon          ","lat          ","time         ","prec         ",&
      "lon          ","lat          ","time         ","ocis         ",&
      "lon          ","lat          ","time         ","olis         "/),(/nvar,nMaxFile2/))

  INTEGER,PUBLIC :: ncid_pres
  INTEGER,PUBLIC :: ncid_prec
  INTEGER,PUBLIC :: ncid_temp
  INTEGER,PUBLIC :: ncid_rhum
  INTEGER,PUBLIC :: ncid_umes
  INTEGER,PUBLIC :: ncid_wind
  INTEGER,PUBLIC :: ncid_swrd
  INTEGER,PUBLIC :: ncid_lwrd

  integer,PUBLIC, ALLOCATABLE :: var_varid    (:,:)
  LOGICAL ::  ALLOC_INIT=.TRUE.
  LOGICAL ::  DEALLOC_FINALIZE=.TRUE.
  PUBLIC :: InitReadNETCDF
  PUBLIC :: FinalizeReadNETCDF
  PUBLIC :: RunReadNETCDF
CONTAINS

 SUBROUTINE InitReadNETCDF(FILE_NAME,icn_data,nFile)
  use netcdf
  IMPLICIT NONE
  ! This is the name of the data file we will read.
  CHARACTER (len = *), INTENT(IN   ) :: FILE_NAME
  INTEGER            , INTENT(IN   ) :: nFile
  CHARACTER (len = *), INTENT(IN   ) :: icn_data
  ! We will learn about the data file and store results in these
  ! program variables.
  INTEGER :: ndims_in
  INTEGER :: nvars_in
  INTEGER :: ngatts_in 
  INTEGER :: unlimdimid_in
  INTEGER :: ivar
  INTEGER :: ncid
  ! Allocate memory.
  IF(ALLOC_INIT)THEN
      IF(TRIM(icn_data) == 'ERA5')THEN
         ALLOCATE(var       (1:nvar,1:nMaxFile))
         ALLOCATE(var_varid (1:nvar,1:nMaxFile))
      ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
         ALLOCATE(var       (1:nvar,1:nMaxFile2))
         ALLOCATE(var_varid (1:nvar,1:nMaxFile2))
      ELSE
          PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
          STOP       
      END IF
      allocate(var_c     (NLONS, NLATS))
      allocate(var_p     (NLONS, NLATS))
      ALLOC_INIT=.FALSE.
  END IF 
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Open the input file. 
  CALL check( nf90_open  (TRIM(FILE_NAME)    , nf90_nowrite, ncid) )
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  !
  ! There are a number of inquiry functions in netCDF which can be
  ! used to learn about an unknown netCDF file. NF90_INQ tells how many
  ! netCDF variables, dimensions, and global attributes are in the
  ! file; also the dimension id of the unlimited dimension, if there
  ! is one.
  call check( nf90_inquire(ncid, ndims_in, nvars_in, ngatts_in, unlimdimid_in))

  !
  !call check( nf90_get_att(ncid, nf90_global, TRIM(history)    , gbl_history))
  !call check( nf90_get_att(ncid, nf90_global, TRIM(Conventions), gbl_Conventions))

  ! Get the varids of the latitude and longitude coordinate variables.
  IF(TRIM(icn_data) == 'ERA5')THEN
      call check( nf90_inq_varid(ncid, VAR_NAME(1,nFile) , var_varid(1,nFile)) )
      call check( nf90_inq_varid(ncid, VAR_NAME(2,nFile) , var_varid(2,nFile)) )
      call check( nf90_inq_varid(ncid, VAR_NAME(3,nFile) , var_varid(3,nFile)) )
  ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
      call check( nf90_inq_varid(ncid, VAR_NAME2(1,nFile) , var_varid(1,nFile)) )
      call check( nf90_inq_varid(ncid, VAR_NAME2(2,nFile) , var_varid(2,nFile)) )
      call check( nf90_inq_varid(ncid, VAR_NAME2(3,nFile) , var_varid(3,nFile)) )
  ELSE
      PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
      STOP       
  END IF
  

PRINT*,'pkubota ',TRIM(FILE_NAME),TRIM(icn_data)
  ! Define the dimensions. The record dimension is defined to have
  ! unlimited length - it can grow as needed. In this example it is
  ! the time dimension.

  ! Read the timer, latitude and longitude data.
  IF(TRIM(icn_data) == 'ERA5')THEN
      call check( nf90_get_var(ncid, var_varid(1,nFile), lons) )
      call check( nf90_get_var(ncid, var_varid(2,nFile), lats) )
  !    call check( nf90_get_var(ncid, var_varid(3,nFile),timer) ) 
  ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
      call check( nf90_get_var(ncid, var_varid(1,nFile), lons) )
      call check( nf90_get_var(ncid, var_varid(2,nFile), lats) )
  !    call check( nf90_get_var(ncid, var_varid(3,nFile),timer) ) 
  ELSE
      PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
      STOP       
  END IF
  PRINT*,'pkubota0 ',TRIM(FILE_NAME),TRIM(icn_data)

  
  ! Define the coordinate variables. We will only define coordinate
  ! variables for lat and lon.  Ordinarily we would need to provide
  ! an array of dimension IDs for each variable's dimensions, but
  ! since coordinate variables only have one dimension, we can
  ! simply provide the address of that dimension ID (lat_dimid) and
  ! similarly for (lon_dimid).

  ! Assign units attributes to the netCDF variables.
  IF(TRIM(icn_data) == 'ERA5')THEN
      call check( nf90_get_att(ncid    ,var_varid(1,nFile), var(1,nFile)%UNITS        , var(1,nFile)%var_units        ) )
      call check( nf90_get_att(ncid    ,var_varid(1,nFile), var(1,nFile)%long_name    , var(1,nFile)%var_long_name    ) )
      call check( nf90_get_att(ncid    ,var_varid(2,nFile), var(2,nFile)%UNITS        , var(2,nFile)%var_units        ) )
      call check( nf90_get_att(ncid    ,var_varid(2,nFile), var(2,nFile)%long_name    , var(2,nFile)%var_long_name    ) )
      call check( nf90_get_att(ncid    ,var_varid(3,nFile), var(3,nFile)%UNITS        , var(3,nFile)%var_units        ) )
      call check( nf90_get_att(ncid    ,var_varid(3,nFile), var(3,nFile)%long_name    , var(3,nFile)%var_long_name    ) )
  ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
      call check( nf90_get_att(ncid    ,var_varid(1,nFile), var(1,nFile)%UNITS        , var(1,nFile)%var_units        ) )
      call check( nf90_get_att(ncid    ,var_varid(1,nFile), var(1,nFile)%long_name    , var(1,nFile)%var_long_name    ) )
      call check( nf90_get_att(ncid    ,var_varid(2,nFile), var(2,nFile)%UNITS        , var(2,nFile)%var_units        ) )
      call check( nf90_get_att(ncid    ,var_varid(2,nFile), var(2,nFile)%long_name    , var(2,nFile)%var_long_name    ) )
      call check( nf90_get_att(ncid    ,var_varid(3,nFile), var(3,nFile)%UNITS        , var(3,nFile)%var_units        ) )
      call check( nf90_get_att(ncid    ,var_varid(3,nFile), var(3,nFile)%long_name    , var(3,nFile)%var_long_name    ) )
  ELSE
      PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
      STOP       
  END IF


PRINT*,'pkubota1 ',TRIM(FILE_NAME),TRIM(icn_data)
  ! The dimids array is used to pass the dimids of the dimensions of
  ! the netCDF variables. Both of the netCDF variables we are creating
  ! share the same four dimensions. In Fortran, the unlimited
  ! dimension must come last on the list of dimids.
  dimids = (/ lon_dimid, lat_dimid, tim_dimid /)

  DO ivar=4,nvar
      !
      ! Get the varids of the pressure and temperature netCDF variables.
      !
      IF(TRIM(icn_data) == 'ERA5')THEN
         call check( nf90_inq_varid(ncid, VAR_NAME(ivar,nFile)  , var_varid(ivar,nFile)) )
      ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
         call check( nf90_inq_varid(ncid, VAR_NAME2(ivar,nFile)  , var_varid(ivar,nFile)) )
      ELSE
         PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
         STOP       
      END IF
      !
      ! Each of the netCDF variables has a "units" attribute. Let's read
      ! them and check them.
      !
!      call check( nf90_get_att(ncid, var_varid(ivar,nFile), var(ivar,nFile)%UNITS        , var(ivar,nFile)%var_units       ) )
!      call check( nf90_get_att(ncid, var_varid(ivar,nFile), var(ivar,nFile)%FillValue    , var(ivar,nFile)%var_FillValue   ) )
      !call check( nf90_get_att(ncid, var_varid(ivar,nFile), var(ivar,nFile)%add_offset   , var(ivar,nFile)%var_add_offset  ) )
      !call check( nf90_get_att(ncid, var_varid(ivar,nFile), var(ivar,nFile)%scale_factor , var(ivar,nFile)%var_scale_factor) )
!      call check( nf90_get_att(ncid, var_varid(ivar,nFile), var(ivar,nFile)%long_name    , var(ivar,nFile)%var_long_name   ) )
  END DO
PRINT*,'pkubota2 ',TRIM(FILE_NAME),TRIM(icn_data)

  IF(TRIM(icn_data) == 'ERA5')THEN
     IF     (nFile ==1) THEN
        ncid_prec=ncid
     ELSE IF(nFile ==2)THEN
        ncid_temp=ncid
     ELSE IF(nFile ==3) THEN
        ncid_rhum=ncid
     ELSE IF(nFile ==4)THEN
        ncid_wind=ncid
     ELSE IF(nFile ==5)THEN
        ncid_swrd=ncid 
     ELSE IF(nFile ==6)THEN
        ncid_lwrd=ncid
     END IF

  ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
     IF     (nFile ==1) THEN
        ncid_pres=ncid
     ELSE IF(nFile ==2)THEN
        ncid_temp=ncid
     ELSE IF(nFile ==3) THEN
        ncid_umes=ncid
     ELSE IF(nFile ==4)THEN
        ncid_wind=ncid
     ELSE IF(nFile ==5)THEN
        ncid_prec=ncid
     ELSE IF(nFile ==6)THEN
        ncid_swrd=ncid
     ELSE IF(nFile ==7)THEN
        ncid_lwrd=ncid
     END IF  
  ELSE
     PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
     STOP       
  END IF


 END SUBROUTINE InitReadNETCDF

 SUBROUTINE RunReadNETCDF(rec,ncid,latic,latip,longc,longp,var_stepc,var_stepp,nFile)
  use netcdf
  IMPLICIT NONE
  INTEGER            , INTENT(IN   ) :: rec
  INTEGER            , INTENT(IN   ) :: nFile
  INTEGER            , INTENT(IN   ) :: ncid
  REAL(KIND=4)       , INTENT(OUT  ) :: latic( NLONS, NLATS)
  REAL(KIND=4)       , INTENT(OUT  ) :: latip( NLONS, NLATS)
  REAL(KIND=4)       , INTENT(OUT  ) :: longc( NLONS, NLATS)
  REAL(KIND=4)       , INTENT(OUT  ) :: longp( NLONS, NLATS)
  REAL(KIND=8)       , INTENT(OUT  ) :: var_stepc( NLONS, NLATS)
  REAL(KIND=8)       , INTENT(OUT  ) :: var_stepp( NLONS, NLATS)
  
  INTEGER  ::ivar,i,j
  !
  ! Read 1 record of NLONS*NLATS values, starting at the beginning 
  ! of the record (the (1, 1, rec) element in the netCDF file).
  !
  count     = (/ NLONS, NLATS, 1 /)
  start     = (/ 1, 1, 1 /)
  start2    = (/ 1, 1, 1 /)
  start_out = (/ 1, 1, 1 /)
  !
  ! Read the surface pressure and temperature data from the file, one
  ! record at a time.
  !
     start(3) = rec
     IF(rec==NTIME)THEN
       start2(3) = rec 
     ELSE 
       start2(3) = rec+1
     END IF
     DO j=1,NLATS
        DO i=1,NLONS
           longp(i,j) =lons(i)
           latip(i,j) =lats(j)
           
           longc(i,j) = longp(i,j)
           latic(i,j) = latip(i,j)
        END DO
     END DO
     DO ivar=4,nvar
        call check( nf90_get_var(ncid, var_varid(ivar,nFile), var_c(1:NLONS, 1:NLATS), start = start  , &
                                 count = count) )
        call check( nf90_get_var(ncid, var_varid(ivar,nFile), var_p(1:NLONS, 1:NLATS), start = start2, &
                                 count = count) )
        DO j=1,NLATS
           DO i=1,NLONS
              var_stepc(i,j) =var_c(i,j)
              var_stepp(i,j) =var_p(i,j)
           END DO
        END DO
     END DO

 END SUBROUTINE RunReadNETCDF

 SUBROUTINE FinalizeReadNETCDF(ncid,nFile)
  use netcdf
  IMPLICIT NONE
  INTEGER            , INTENT(IN   ) :: nFile
  INTEGER            , INTENT(IN   ) :: ncid
  IF(DEALLOC_FINALIZE)THEN
     DEALLOCATE(var)
     DEALLOCATE(var_varid)
     DEALLOCATE(var_c  )
     DEALLOCATE(var_p )
     DEALLOC_FINALIZE=.FALSE.
  END IF 
  ! Close the file. This frees up any internal netCDF resources
  ! associated with the file.
  IF(nFile==1) call check( nf90_close(ncid) )
 END SUBROUTINE FinalizeReadNETCDF


  subroutine check(status)
  use netcdf
    integer, intent ( in) :: status    
    if(status /= nf90_noerr) then 
      stop 2
    end if
  end subroutine check  
END MODULE ReadNETCDF
