! 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

! : pres_temp_4D_rd.f90,v 1.8 2008/02/19 16:32:21 ed Exp $

program pres_temp_4D_rd
  use netcdf
  implicit none
  
  ! This is the name of the data file we will read.
  character (len = *), parameter :: FILE_NAME = "ta.mes28.nc"
  character (len = *), parameter :: FILE_NAME1 = "ws.mes28.nc"
  integer :: ncid

  ! We are reading 4D data, a 12 x 6 x 2 lon-lat-lvl grid, with 2
  ! timesteps of data.
  integer, parameter :: NDIMS = 4, NRECS = 112
  integer, parameter :: NLVLS = 23, NLATS = 160, NLONS = 320
  character (len = *), parameter :: LVL_NAME = "lev"
  character (len = *), parameter :: LAT_NAME = "lat"
  character (len = *), parameter :: LON_NAME = "lon"
  character (len = *), parameter :: REC_NAME = "time"
  integer :: lvl_dimid, lon_dimid, lat_dimid, rec_dimid

  ! The start and count arrays will tell the netCDF library where to
  ! read our data.
  integer :: start(NDIMS), count(NDIMS)

  ! 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 :: lvls(NLVLS), lats(NLATS), lons(NLONS), tims(NRECS)
  integer :: lon_varid, lat_varid, lvl_varid, tim_varid

  ! We will read surface temperature and pressure fields. In netCDF
  ! terminology these are called "variables."
  character (len = *), parameter :: TEMP_NAME="ta"
  integer :: temp_varid
  integer :: dimids(NDIMS)

  ! We recommend that each variable carry a "units" attribute.
  character (len = *), parameter :: UNITS = "units"
  character (len = *), parameter :: LVL_UNITS = "hPa", TEMP_UNITS = "celsius"
  character (len = *), parameter :: LAT_UNITS = "degrees_north"
  character (len = *), parameter :: LON_UNITS = "degrees_east"
  character*100 TIM_UNITS

  ! Program variables to hold the data we will read in. We will only
  ! need enough space to hold one timestep of data; one record.
  real :: temp_in(NLONS, NLATS, NLVLS),temp_out(NLONS, NLATS, NLVLS,NRECS)

  ! Loop indices
  integer :: lvl, lat, lon, rec, i

  ! To check the units attributes.
  character*80 temp_units_in
  character*80 lat_units_in, lon_units_in

  ! Open the file. 
  call check( nf90_open(FILE_NAME, nf90_nowrite, ncid) )

  ! Get the varids of the latitude and longitude coordinate variables.
  call check( nf90_inq_varid(ncid, LVL_NAME, lvl_varid) )
  call check( nf90_inq_varid(ncid, LAT_NAME, lat_varid) )
  call check( nf90_inq_varid(ncid, LON_NAME, lon_varid) )
  call check( nf90_inq_varid(ncid, REC_NAME, tim_varid) )

  ! Read the latitude and longitude data.
  call check( nf90_get_var(ncid, lvl_varid, lvls) )
  call check( nf90_get_var(ncid, lat_varid, lats) )
  call check( nf90_get_var(ncid, lon_varid, lons) )
  call check( nf90_get_var(ncid, tim_varid, tims) )

  ! Assign units attributes to coordinate variables.
  call check( nf90_get_att(ncid, tim_varid, UNITS, TIM_UNITS) )

  ! Get the varids of the pressure and temperature netCDF variables.
  call check( nf90_inq_varid(ncid, TEMP_NAME, temp_varid) )

  ! Read 1 record of NLONS*NLATS*NLVLS values, starting at the beginning 
  ! of the record (the (1, 1, 1, rec) element in the netCDF file).
  count = (/ NLONS, NLATS, NLVLS, 1 /)
  start = (/ 1, 1, 1, 1 /)

  ! Read the surface pressure and temperature data from the file, one
  ! record at a time.
  do rec = 1, NRECS
     start(4) = rec
     call check( nf90_get_var(ncid, temp_varid, temp_in, start, count) )
     
     i = 0
     do lvl = 1, NLVLS
!     print*,lvls(lvl)
        do lat = 1, NLATS
           do lon = 1, NLONS

! a variavel temp_in(lon, lat, lvl) deve ser manipulada aqui. Ela e a variavel de input
! a variavel temp_out(lon, lat,lvl, rec) e a variavel final a ser calculada e escrita no novo arquivo. 

! PRESSAO DE SATURACAO DO VAPOR (es) i
            if(temp_in(lon,lat,lvl)==-999)then
            temp_out(lon,lat,lvl,rec)=temp_in(lon,lat,lvl)
            else
            
              if(temp_in(lon,lat,lvl)>=273.15)then
                temp_out(lon,lat,lvl,rec)=6.1078*(10**((7.5*(temp_in(lon,lat,lvl)-273.15))/(237.3+(temp_in(lon,lat,lvl)-273.15))))
              else
                temp_out(lon,lat,lvl,rec)=6.1078*(10**((9.5*(temp_in(lon,lat,lvl)-273.15))/(265.5+(temp_in(lon,lat,lvl)-273.15))))
              endif
! RAZAO DE MISTURA DE SATURACAO (Ws)
              temp_out(lon,lat,lvl,rec)=0.622*(temp_out(lon,lat,lvl,rec)/(lvls(lvl)-temp_out(lon,lat,lvl,rec)))
! DIVIDINDO POR 100 PRA DEPOIS FAZER UR=W/Ws. COMO JA ESTOU DIVIDINDO POR 100 AQUI, DEPOIS NAO PRECISO MULTIPLICAR UR POR 100
              temp_out(lon,lat,lvl,rec)=temp_out(lon,lat,lvl,rec)/100
              i = i + 1
            endif
           end do
        end do
     end do
     ! next record
  end do

         
  ! Close the file. This frees up any internal netCDF resources
  ! associated with the file.
  call check( nf90_close(ncid) )

  ! If we got this far, everything worked as expected. Yipee! 
  print *,"*** SUCCESS reading example file ", FILE_NAME, "!"


! PARTE QUE CRIA O OUTRO ARQUIVO

  ! Create the file. 
  call check( nf90_create(FILE_NAME1, nf90_clobber, ncid) )

  ! 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.
  call check( nf90_def_dim(ncid, LVL_NAME, NLVLS, lvl_dimid) )
  call check( nf90_def_dim(ncid, LAT_NAME, NLATS, lat_dimid) )
  call check( nf90_def_dim(ncid, LON_NAME, NLONS, lon_dimid) )
!  call check( nf90_def_dim(ncid, REC_NAME, NRECS, rec_dimid) )
  call check( nf90_def_dim(ncid, REC_NAME, NF90_UNLIMITED, rec_dimid) )

  ! 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).
  call check( nf90_def_var(ncid, LVL_NAME, NF90_REAL, lvl_dimid, lvl_varid) )
  call check( nf90_def_var(ncid, LAT_NAME, NF90_REAL, lat_dimid, lat_varid) )
  call check( nf90_def_var(ncid, LON_NAME, NF90_REAL, lon_dimid, lon_varid) )
  call check( nf90_def_var(ncid, REC_NAME, NF90_REAL, rec_dimid, tim_varid) )

  ! Assign units attributes to coordinate variables.
  call check( nf90_put_att(ncid, lvl_varid, UNITS, LVL_UNITS) )
  call check( nf90_put_att(ncid, lat_varid, UNITS, LAT_UNITS) )
  call check( nf90_put_att(ncid, lon_varid, UNITS, LON_UNITS) )
  call check( nf90_put_att(ncid, tim_varid, UNITS, TIM_UNITS) )

  ! 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, lvl_dimid, rec_dimid /)

  ! Define the netCDF variables for the pressure and temperature data.
  call check( nf90_def_var(ncid, TEMP_NAME, NF90_REAL, dimids, temp_varid) )

  ! Assign units attributes to the netCDF variables.
  call check( nf90_put_att(ncid, temp_varid, UNITS, TEMP_UNITS) )
  call check( nf90_put_att(ncid, temp_varid, '_FillValue', -999.) )

  ! End define mode.
  call check( nf90_enddef(ncid) )

  ! Write the coordinate variable data. This will put the latitudes
  ! and longitudes of our data grid into the netCDF file.
  call check( nf90_put_var(ncid, lvl_varid, lvls) )
  call check( nf90_put_var(ncid, lat_varid, lats) )
  call check( nf90_put_var(ncid, lon_varid, lons) )
  call check( nf90_put_var(ncid, tim_varid, tims) )

  ! These settings tell netcdf to write one timestep of data. (The
  ! setting of start(4) inside the loop below tells netCDF which
  ! timestep to write.)
  count = (/ NLONS, NLATS, NLVLS, NRECS /)
!  count = (/ NLONS, NLATS, NLVLS, 1 /)
  start = (/ 1, 1, 1, 1 /)

  ! Write the pretend data. This will write our surface pressure and
  ! surface temperature data. The arrays only hold one timestep worth
  ! of data. We will just rewrite the same data for each timestep. In
  ! a real :: application, the data would change between timesteps.
  do rec = 1, 1
     start(4) = rec
     call check( nf90_put_var(ncid, temp_varid, temp_out, start, count) )
!     call check( nf90_put_var(ncid, temp_varid, temp_out, start = start, &
!                              count = count) )
  end do

  ! Close the file. This causes netCDF to flush all buffers and make
  ! sure your data are really written to disk.
  call check( nf90_close(ncid) )

  print *,"*** SUCCESS writing example file ", FILE_NAME1, "!"


contains
  subroutine check(status)
    integer, intent ( in) :: status
    
    if(status /= nf90_noerr) then 
      print *, trim(nf90_strerror(status))
      stop 2
    end if
  end subroutine check  
end program pres_temp_4D_rd