MODULE MOD_NETCDF_xice
  use netcdf

 IMPLICIT NONE
   PRIVATE
   INTEGER, PARAMETER :: &
            r4 = SELECTED_REAL_KIND(6)  ! Kind for 32-bits Real Numbers
   INTEGER, PARAMETER :: &
            r8 = SELECTED_REAL_KIND(15) ! Kind for 64-bits Real Numbers
   INTEGER, PARAMETER    :: nFile      =   1
   integer, parameter    :: nMaxFile   =   1
   integer, parameter    :: nvar       =   4
   integer, parameter    :: natrb      =   8
   INTEGER, PARAMETER    :: MAX_ATT_LEN= 800
   INTEGER, PUBLIC       :: xsize_xice, ysize_xice , tsize_xice
   INTEGER               :: xsize, ysize , tsize

   INTEGER, PARAMETER, PUBLIC :: nVarNetcdf_xice(1:nvar)=RESHAPE((/&
       1                    ,2                     ,3                     ,4                     /),(/nvar/))

  CHARACTER (len = 20), PARAMETER,PUBLIC :: VAR_NAME(1:nvar,1:nMaxFile)=RESHAPE((/&
      "time                ","lat                 ","lon                 ","icec                "/),(/nvar,nMaxFile/))

  ! We recommend that each variable carry a "units" attribute.
  TYPE :: attribute_latitude
     character (len = 13) :: long_name               = "long_name" ;
     character (len = 13) :: standard_name           = "standard_name" ;
     character (len = 13) :: units                   = "units" ;
     character (len = 13) :: actual_range            = "actual_range" ;
     character (len = 13) :: axis                    = "axis" ;

     character (len = 15) :: var_long_name            != "Latitude" ;
     character (len = 99) :: var_standard_name        != "latitude" ;
     character (len = 15) :: var_units                != "degrees_north" ;
     REAL(KIND=r4)        :: var_actual_range(2)      != "-89.875f, 89.875f " ;
     character (len = 99) :: var_axis                 != "y" ;
  END TYPE attribute_latitude
  TYPE(attribute_latitude) :: Coord_lat


  ! We recommend that each variable carry a "units" attribute.
  TYPE :: attribute_Longitude
     character (len = 13) :: long_name               = "long_name" ;
     character (len = 13) :: standard_name           = "standard_name" ;
     character (len = 13) :: units                   = "units" ;
     character (len = 13) :: actual_range            = "actual_range" ;
     character (len = 13) :: axis                    = "axis" ;

     character (len = 15) :: var_long_name            != "Longitude" ;
     character (len = 99) :: var_standard_name        != "longitude" ;
     character (len = 15) :: var_units                != "degrees_east" ;
     REAL(KIND=r4)        :: var_actual_range(2)      != " 0.125f, 359.875f " ;
     character (len = 99) :: var_axis                 != "X" ;

  END TYPE attribute_Longitude
  TYPE(attribute_Longitude) :: Coord_lon
 
  ! We recommend that each variable carry a "units" attribute.
  TYPE :: attribute_time
     character (len = 13) :: long_name             = "long_name"
     character (len = 13) :: units                 = "units" ;
     character (len = 13) :: delta_t               = "delta_t"
     character (len = 13) :: avg_period            = "avg_period"  
     character (len = 13) :: axis                  = "axis" ;
     character (len = 13) :: actual_range          = "actual_range" ;

     character (len = 99) :: var_long_name           != "Time" ;
     character (len = 99) :: var_units               != "days since 1800-01-01 00:00:00" ;
     character (len = 99) :: var_delta_t             != "0000-00-01 00:00:00" ;
     character (len = 99) :: var_avg_period          != "0000-00-01 00:00:00"    
     character (len = 99) :: var_axis                != "T" ;
     REAL(KIND=r8)        :: var_actual_range(2)     != "81449., 81506." ;
  END TYPE attribute_time
  TYPE(attribute_time) :: Coordtime


  ! We recommend that each variable carry a "units" attribute.
  TYPE :: attribute_sfc_xice
     character (len = 13) :: long_name        = "long_name" ;
     character (len = 13) :: units            = "units" ;
     character (len = 13) :: valid_range      = "valid_range" ;
     character (len = 13) :: missing_value    = "missing_value" ;
     character (len = 13) :: precision        = "precision"
     character (len = 13) :: dataset          = "dataset"
     character (len = 13) :: var_desc         = "var_desc"
     character (len = 13) :: level_desc       = "level_desc"
     character (len = 13) :: statistic        = "statistic"
     character (len = 13) :: parent_stat      = "parent_stat"
     character (len = 13) :: actual_range     = "actual_range"

     character (len = 99)  :: var_long_name       != "Daily Sea Surface Temperature" ;
     character (len = 99)  :: var_units           != "degC" ;
     REAL(KIND=r4)         :: var_valid_range(2)  != -3.f, 45.f  ;
     REAL(KIND=r4)         :: var_missing_value   != -9.96921e+36f ;
     REAL(KIND=r4)         :: var_precision       ! = 2.f
     character (len = 99)  :: var_dataset         ! = "NOAA High-resolution Blended Analysis"
     character (len = 99)  :: var_var_desc        ! = "Sea Surface Temperature"
     character (len = 99)  :: var_level_desc      ! = "Surface"
     character (len = 99)  :: var_statistic       ! = "Mean"
     character (len = 99)  :: var_parent_stat     ! = "Individual Observations"
     REAL(KIND=r4)         :: var_actual_range(2) ! =-1.8f, 33.45f ;
   END TYPE attribute_sfc_xice
  TYPE(attribute_sfc_xice),PUBLIC :: sfc_xice


   CHARACTER (len = 526) :: FILE_NAME
   integer               :: ncid
   INTEGER               :: ndims_in  
   INTEGER               :: nvars_in
   INTEGER               :: ngatts_in 
   INTEGER               :: unlimdimid_in
   CHARACTER(LEN=50)     :: xname, yname,zname,zname2, tname
   INTEGER               :: ivar 
   ! We will read surface temperature and pressure fields. In netCDF
   ! terminology these are called "variables."
!   CHARACTER (len = 26), PARAMETER :: ATR_NAME(1:natrb)=RESHAPE((/&
!        "Conventions   ",&
!        "title         ",&
!        "institution   ",&
!        "source        ",&
!        "history       ",&
!        "dataset_title ",&
!        "References    ",&
!        "comment       "/),(/natrb/))

   ! We will read surface temperature and pressure fields. In netCDF
   ! terminology these are called "variables."
   CHARACTER (len = 26), PARAMETER :: ATR_NAME(1:natrb)=RESHAPE((/&
        "Conventions   ",&
        "title         ",&
        "institution   ",&
        "source        ",&
        "References    ",&
        "dataset_title ",&
        "version       ",&
        "comment       "/),(/natrb/))

  TYPE :: GBL_attribute
    INTEGER                       :: IntegerGlobal_tm       !"= "UNCLASSIFIED" ;
    CHARACTER (len = MAX_ATT_LEN) :: StringGlobal(natrb)    !"= "UNCLASSIFIED" ;
    INTEGER                       :: IntegerGlobal_im       !"= "UNCLASSIFIED" ;
    INTEGER                       :: IntegerGlobal_jm       !"= "UNCLASSIFIED" ;
  END TYPE GBL_attribute
  TYPE(GBL_attribute), ALLOCATABLE :: Agbl(:)


  integer, ALLOCATABLE :: var_varid    (:,:)
  integer, ALLOCATABLE :: start        (:)
  INTEGER, ALLOCATABLE :: count        (:)
  REAL (KIND=r4),            ALLOCATABLE  :: buff(:)
  REAL (KIND=r8),PUBLIC    , ALLOCATABLE  :: lon0_xice(:)
  REAL (KIND=r8),PUBLIC    , ALLOCATABLE  :: lat0_xice(:)
  
  PUBLIC :: Init_Netcdf_Parameter_xice 
  PUBLIC :: READ_Netcdf_Field_xice
  PUBLIC :: READ_Netcdf_Field_xice2
  PUBLIC :: Finalize_Netcdf_Parameter_xice
  
CONTAINS

 SUBROUTINE Init_Netcdf_Parameter_xice(FILE_NAME)
   IMPLICIT NONE
   CHARACTER(LEN=*), INTENT(IN   ) :: FILE_NAME
   INTEGER :: fieldsize,im,jm,tm,ij,i,j
   PRINT*,TRIM(FILE_NAME)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Allocate memory.
   !  IF(ALLOC_INIT)THEN
   IF (.not. Allocated( Agbl      ) )   ALLOCATE(Agbl      (1:nMaxFile))
   If (.not. Allocated( var_varid ) )   ALLOCATE(var_varid (1:nvar,1:nMaxFile))
   !     ALLOC_INIT=.FALSE.
   !  END IF 
   ! Open the input file. 
   CALL check( nf90_open  (TRIM(FILE_NAME)    , nf90_nowrite, ncid) )

   PRINT*,'passei 0'

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !
   !
   ! There are a number of inquiry functions in netCDF which can be
   ! used to learn about an unknown netCDF file. NF90_INQ tells how maysize
   ! 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))


   PRINT*,'passei 1'
   
   !Inquire about the dimensions
   !:-------:-------:-------:-------:-------:-------:-------:
   CALL check(nf90_inquire_dimension(ncid,1,tname ,tsize))
   CALL check(nf90_inquire_dimension(ncid,2,yname ,ysize))
   CALL check(nf90_inquire_dimension(ncid,3,xname ,xsize))
   PRINT*,'passei 2'
   xsize_xice =xsize
   ysize_xice =ysize
   tsize_xice =tsize
   Agbl(1)%IntegerGlobal_im=xsize
   Agbl(1)%IntegerGlobal_jm=ysize
   Agbl(1)%IntegerGlobal_tm=tsize
   !
   ! Get the varids of the pressure and temperature netCDF variables.
   !
   !:-------:-------:-------:-------:-------:-------:-------:
   call check( nf90_get_att(ncid, nf90_global,  ATR_NAME(1)    , Agbl(nFile)%StringGlobal(1)))
   call check( nf90_get_att(ncid, nf90_global,  ATR_NAME(2)    , Agbl(nFile)%StringGlobal(2)))
   call check( nf90_get_att(ncid, nf90_global,  ATR_NAME(3)    , Agbl(nFile)%StringGlobal(3)))
   call check( nf90_get_att(ncid, nf90_global,  ATR_NAME(4)    , Agbl(nFile)%StringGlobal(4)))
   call check( nf90_get_att(ncid, nf90_global,  ATR_NAME(5)    , Agbl(nFile)%StringGlobal(5)))
   call check( nf90_get_att(ncid, nf90_global,  ATR_NAME(6)    , Agbl(nFile)%StringGlobal(6)))
   call check( nf90_get_att(ncid, nf90_global,  ATR_NAME(7)    , Agbl(nFile)%StringGlobal(7)))
   call check( nf90_get_att(ncid, nf90_global,  ATR_NAME(8)    , Agbl(nFile)%StringGlobal(8)))

   PRINT*,'passei 3'

   ! Get the varids of the latitude and longitude coordinate variables.

   DO ivar=1,nvar
      call check( nf90_inq_varid(ncid, VAR_NAME(ivar,nFile) , var_varid(ivar,nFile)) )
   END DO 

   PRINT*,'passei 4'

   ! 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.

   ! We recommend that each variable carry a "units" attribute.

   call check( nf90_get_att(ncid    ,var_varid(2,nFile),  Coord_lat%long_name        , Coord_lat%var_long_name     ) )
   call check( nf90_get_att(ncid    ,var_varid(2,nFile),  Coord_lat%standard_name    , Coord_lat%var_standard_name ) )
   call check( nf90_get_att(ncid    ,var_varid(2,nFile),  Coord_lat%units            , Coord_lat%var_units         ) )
   call check( nf90_get_att(ncid    ,var_varid(2,nFile),  Coord_lat%actual_range     , Coord_lat%var_actual_range  ) )
   call check( nf90_get_att(ncid    ,var_varid(2,nFile),  Coord_lat%axis             , Coord_lat%var_axis          ) )

   PRINT*,'passei 5'
   call check( nf90_get_att(ncid    ,var_varid(3,nFile),  Coord_lon%long_name        , Coord_lon%var_long_name     ) )
   call check( nf90_get_att(ncid    ,var_varid(3,nFile),  Coord_lon%standard_name    , Coord_lon%var_standard_name ) )
   call check( nf90_get_att(ncid    ,var_varid(3,nFile),  Coord_lon%units            , Coord_lon%var_units         ) )
   call check( nf90_get_att(ncid    ,var_varid(3,nFile),  Coord_lon%actual_range     , Coord_lon%var_actual_range  ) )
   call check( nf90_get_att(ncid    ,var_varid(3,nFile),  Coord_lon%axis             , Coord_lon%var_axis          ) )

   PRINT*,'passei 6'

   call check( nf90_get_att(ncid    ,var_varid(1,nFile),  Coordtime%long_name      , Coordtime%var_long_name         ) )
   call check( nf90_get_att(ncid    ,var_varid(1,nFile),  Coordtime%units          , Coordtime%var_units             ) )
   call check( nf90_get_att(ncid    ,var_varid(1,nFile),  Coordtime%delta_t        , Coordtime%var_delta_t           ) )
   call check( nf90_get_att(ncid    ,var_varid(1,nFile),  Coordtime%avg_period     , Coordtime%var_avg_period        ) )
   call check( nf90_get_att(ncid    ,var_varid(1,nFile),  Coordtime%axis           , Coordtime%var_axis              ) )
   call check( nf90_get_att(ncid    ,var_varid(1,nFile),  Coordtime%actual_range   , Coordtime%var_actual_range      ) )
   PRINT*,'passei 7'
     
   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%long_name          ,  sfc_xice%var_long_name      ))
   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%units              ,  sfc_xice%var_units          ))
   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%valid_range        ,  sfc_xice%var_valid_range    ))
   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%missing_value      ,  sfc_xice%var_missing_value  ))
   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%precision          ,  sfc_xice%var_precision      ))
   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%dataset            ,  sfc_xice%var_dataset        ))
   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%var_desc           ,  sfc_xice%var_var_desc       ))
   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%level_desc         ,  sfc_xice%var_level_desc     ))
   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%statistic          ,  sfc_xice%var_statistic      ))

   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%parent_stat        ,  sfc_xice%var_parent_stat    ))
   call check( nf90_get_att(ncid    ,var_varid(4,nFile),  sfc_xice%actual_range       ,  sfc_xice%var_actual_range   ))

   !PRINT*,'Longitude 8'
   im=Agbl(1)%IntegerGlobal_im
   jm=Agbl(1)%IntegerGlobal_jm
   tm=Agbl(1)%IntegerGlobal_tm
   fieldsize = im
   If (.not. Allocated( count      ) )   ALLOCATE(count      (3))
   If (.not. Allocated( start      ) )   ALLOCATE(start      (3))
   If (.not. Allocated( buff       ) )   ALLOCATE(buff(fieldsize))
   If (.not. Allocated( lon0_xice   ) )   ALLOCATE(lon0_xice(fieldsize))

   ! Fields DIMS
   count     = (/    im,     1,   1 /)
   start     = (/     1,     1,   1 /)
   start(3)  = 1
   call check2( nf90_get_var(ncid, var_varid(3,nFile),buff,start = start  , count = count))
   WRITE (UNIT=*, FMT='(A20,1P2G12.5)') TRIM(VAR_NAME(3,nFile)),MINVAL(buff),MAXVAL(buff)

   ij=0
   DO i=1,im
      ij=ij+1
      lon0_xice(i) = REAL(buff(ij),KIND=r8)
   END DO
   If (Allocated( count      ) )   DEALLOCATE(count      )
   If (Allocated( start      ) )   DEALLOCATE(start      )
   If (Allocated( buff       ) )   DEALLOCATE(buff      )


   PRINT*,'latitude 9'
   fieldsize = jm
   If (.not. Allocated( count      ) )   ALLOCATE(count      (3))
   If (.not. Allocated( start      ) )   ALLOCATE(start      (3))
   If (.not. Allocated( buff       ) )   ALLOCATE(buff(fieldsize))
   If (.not. Allocated( lat0_xice   ) )   ALLOCATE(lat0_xice(fieldsize))

   ! Fields DIMS
   count     = (/    jm,     1,   1 /)
   start     = (/     1,     1,   1 /)
   start(3)  = 1
   call check2( nf90_get_var(ncid, var_varid(2,nFile),buff,start = start  , count = count))
   WRITE (UNIT=*, FMT='(A20,1P2G12.5)') TRIM(VAR_NAME(2,nFile)),MINVAL(buff),MAXVAL(buff)
   ij=0
   DO i=1,jm
      ij=ij+1
      lat0_xice(i) = REAL(buff(ij),KIND=r8)
   END DO
   If (Allocated( count      ) )   DEALLOCATE(count      )
   If (Allocated( start      ) )   DEALLOCATE(start      )
   If (Allocated( buff       ) )   DEALLOCATE(buff       )

 END SUBROUTINE Init_Netcdf_Parameter_xice

 SUBROUTINE READ_Netcdf_Field_xice2(SoilGdas,irec_time,nfield)
  IMPLICIT NONE
  REAL(KIND=r8), INTENT(OUT  ) :: SoilGdas(:,:)
  INTEGER      , INTENT(IN   ) :: nfield 
  INTEGER      , INTENT(IN   ) :: irec_time

  INTEGER :: fieldsize,im,jm,tm,i,j,ij
  im=Agbl(1)%IntegerGlobal_im
  jm=Agbl(1)%IntegerGlobal_jm
  tm=Agbl(1)%IntegerGlobal_tm
  fieldsize = im*jm
  If (.not. Allocated( count      ) )   ALLOCATE(count      (3))
  If (.not. Allocated( start      ) )   ALLOCATE(start      (3))
  If (.not. Allocated( buff       ) )   ALLOCATE(buff(fieldsize))

  ! Fields DIMS
  count     = (/    im,    jm,   1 /)
  start     = (/     1,     1,   1 /)
  start(3) = irec_time
  call check2( nf90_get_var(ncid, var_varid(nfield,nFile),buff,start = start  , count = count))
  WRITE (UNIT=*, FMT='(A20,I5,1P3G12.5)') TRIM(VAR_NAME(nfield,1)),start(3), buff(1),MINVAL(buff),MAXVAL(buff)
  ij=0
  DO j=1,jm
     DO i=1,im
       ij=ij+1
       SoilGdas(i,j) = REAL(buff(ij),KIND=r8)
     END DO
  END DO
  If (Allocated( count      ) )   DEALLOCATE(count      )
  If (Allocated( start      ) )   DEALLOCATE(start      )
 END SUBROUTINE READ_Netcdf_Field_xice2 
 
 SUBROUTINE READ_Netcdf_Field_xice(SoilGdas,nfield)

  IMPLICIT NONE
  REAL(KIND=r8), INTENT(OUT  ) :: SoilGdas(:,:)
  INTEGER      , INTENT(IN   ) :: nfield 
  INTEGER :: fieldsize,im,jm,tm,i,j,ij
  im=Agbl(1)%IntegerGlobal_im
  jm=Agbl(1)%IntegerGlobal_jm
  tm=Agbl(1)%IntegerGlobal_tm
  fieldsize = im*jm
  If (.not. Allocated( count      ) )   ALLOCATE(count      (3))
  If (.not. Allocated( start      ) )   ALLOCATE(start      (3))
  If (.not. Allocated( buff       ) )   ALLOCATE(buff(fieldsize))

  ! Fields DIMS
  count     = (/    im,    jm,   1 /)
  start     = (/     1,     1,   1 /)
  start(3) = 1
  call check( nf90_get_var(ncid, var_varid(nfield,nFile),buff,start = start  , count = count))
  !PK WRITE (UNIT=*, FMT='(A20,I5,1P3G12.5)') TRIM(VAR_NAME(nfield,1)),start(3), buff(1),MINVAL(buff),MAXVAL(buff)
  ij=0
  DO j=1,jm
     DO i=1,im
       ij=ij+1
       SoilGdas(i,j) = REAL(buff(ij),KIND=r8)
     END DO
  END DO
  If (Allocated( count      ) )   DEALLOCATE(count      )
  If (Allocated( start      ) )   DEALLOCATE(start      )
 END SUBROUTINE READ_Netcdf_Field_xice
 
 
 SUBROUTINE Finalize_Netcdf_Parameter_xice()
   IMPLICIT NONE
  !Close netcdf files
   call check( nf90_close(ncid) )

 END SUBROUTINE Finalize_Netcdf_Parameter_xice

  SUBROUTINE check(status)
  USE netcdf
    INTEGER, INTENT ( in) :: status
    IF(status /= nf90_noerr) THEN 
      STOP 2
    END IF
  END SUBROUTINE check  

  SUBROUTINE check2(status)
  USE netcdf
    INTEGER, INTENT ( in) :: status
    IF(status /= nf90_noerr) THEN 
      STOP 'no data'
    END IF
  END SUBROUTINE check2 

END MODULE MOD_NETCDF_xice
