MODULE read_data_PCM_mod

use netcdf, only: nf90_open, NF90_NOWRITE, nf90_noerr, nf90_strerror, &
                  nf90_get_var, nf90_inq_varid, nf90_inq_dimid,       &
                  nf90_inquire_dimension, nf90_close

implicit none

character(13), parameter :: modname = 'read_data_PCM'
character(256)           :: msg      ! for reading
integer                  :: fID, vID ! for reading

!=======================================================================
contains
!=======================================================================

SUBROUTINE read_data_PCM(filename1,filename2,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_PCM,ps_timeseries,ps_avg,tsurf_avg_yr1,tsurf_avg, &
                         tsoil_avg,tsoil_timeseries,min_co2_ice,min_h2o_ice,q_co2,q_h2o,watersurf_density_avg,watersoil_density_timeseries)
use comsoil_h,             only: nsoilmx
use comsoil_h_PEM,         only: soil_pem
use constants_marspem_mod, only: m_co2, m_noco2

implicit none

!=======================================================================
!
! Purpose: Read initial confitions file from the PCM
!
! Authors: JBC
!=======================================================================

include "dimensions.h"

!=======================================================================
! Inputs:
character(*), intent(in) :: filename1, filename2                ! File names
integer,      intent(in) :: timelen                             ! Number of times stored in the file
integer,      intent(in) :: iim_input, jjm_input, ngrid, nslope ! Number of points in the lat x lon dynamical grid, number of subgrid slopes
! Ouputs
real, dimension(ngrid,nslope,2),               intent(out) :: min_co2_ice                  ! Minimum of co2 ice per slope of the year [kg/m^2]
real, dimension(ngrid,nslope,2),               intent(out) :: min_h2o_ice                  ! Minimum of h2o ice per slope of the year [kg/m^2]
real, dimension(ngrid,timelen),                intent(out) :: vmr_co2_PCM                  ! Grid points x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
real, dimension(ngrid,timelen),                intent(out) :: q_co2                        ! CO2 mass mixing ratio in the first layer [kg/m^3]
real, dimension(ngrid,timelen),                intent(out) :: q_h2o                        ! H2O mass mixing ratio in the first layer [kg/m^3]
real, dimension(ngrid,timelen),                intent(out) :: ps_timeseries                ! Surface pressure timeseries [Pa]
real, dimension(ngrid),                        intent(out) :: ps_avg                       ! Averaged surface pressure [K]
real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg                    ! Averaged surface temperature [K]
real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg_yr1                ! Averaged surface temperature for year 1 [K]
real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_avg                    ! Averaged soil temperature for year 2 [K]
real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_timeseries             ! Soil temperature timeseries [K]
real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_avg        ! Water density at the surface [kg/m^3]
real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density_timeseries ! Water density timeseries in the soil layer [kg/m^3]
! Local variables
integer                               :: i, j, l, islope                       ! Loop variables
real                                  :: A, B                                  ! Intermediate variables to compute the mean molar mass of the layer
character(:), allocatable             :: num                                   ! For reading sloped variables
real, dimension(:,:),     allocatable :: var2_read                             ! Variables for reading (2 dimensions)
real, dimension(:,:,:),   allocatable :: var3_read_1, var3_read_2, var3_read_3 ! Variables for reading (3 dimensions)
real, dimension(:,:,:,:), allocatable :: var4_read                             ! Variables for reading (4 dimensions)
!-----------------------------------------------------------------------

!------------------------------- Year 1 --------------------------------
! Open the NetCDF file
write(*,*) "Opening "//filename1//"..."
call error_msg(NF90_OPEN(filename1,NF90_NOWRITE,fID),"open",filename1)

! Download the data from the file
allocate(var2_read(iim_input + 1,jjm_input + 1),var3_read_1(iim_input + 1,jjm_input + 1,timelen),var3_read_2(iim_input + 1,jjm_input + 1,timelen))

if (nslope == 1) then ! No slope
    allocate(character(0) :: num)
else ! We use slopes
    allocate(character(8) :: num)
endif

do islope = 1,nslope
    if (nslope /= 1) then
        write(num,'(i2.2)') islope
        num = '_slope'//num
    endif

    ! CO2 ice
    !--------
    call get_var3("co2ice"//num,var3_read_1)
    where (var3_read_1 < 0.) var3_read_1 = 0.
    write(*,*) "Data for co2_ice"//num//" downloaded."
#ifndef CPP_STD
    call get_var3("perennial_co2ice"//num,var3_read_2)
    write(*,*) "Data for perennial_co2ice"//num//" downloaded."
#endif

    ! Compute the minimum over the year for each point
    var2_read = minval(var3_read_1 + var3_read_2,3)
#ifndef CPP_1D
    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,1))
#else
    min_co2_ice(1,islope,1) = var2_read(1,1)
#endif
    write(*,*) "Min of co2_ice"//num//" computed."

    ! H2O ice
    !--------
    call get_var3("h2o_ice_s"//num,var3_read_1)
    where (var3_read_1 < 0.) var3_read_1 = 0.
    write(*,*) "Data for h2o_ice_s"//num//" downloaded."
#ifndef CPP_STD
    call get_var3("watercap"//num,var3_read_2)
    write(*,*) "Data for watercap"//num//" downloaded."
#endif

    ! Compute the minimum over the year for each point
    var2_read = minval(var3_read_1 + var3_read_2,3)
#ifndef CPP_1D
    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,1))
#else
    min_h2o_ice(1,islope,1) = var2_read(1,1)
#endif
    write(*,*) "Min of h2o_ice"//num//" computed."

    ! Tsurf
    !------
    call get_var3("tsurf"//num,var3_read_1)
    write(*,*) "Data for tsurf"//num//" downloaded."

    ! Compute average over the year for each point
    var2_read = sum(var3_read_1,3)/timelen
#ifndef CPP_1D
    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,islope))
#else
    tsurf_avg_yr1(1,islope) = var2_read(1,1)
#endif
    write(*,*) "Average of tsurf"//num//" computed."
enddo

! Close the NetCDF file
call error_msg(nf90_close(fID),"close",filename1)
write(*,*) '> "'//filename1//'" processed!'

!------------------------------- Year 2 --------------------------------
! Open the NetCDF file
write(*,*) "Opening "//filename2//"..."
call error_msg(NF90_OPEN(filename2,NF90_NOWRITE,fID),"open",filename2)

! Download the data from the file
allocate(var3_read_3(iim_input + 1,jjm_input + 1,nsoilmx),var4_read(iim_input + 1,jjm_input + 1,nsoilmx,timelen))

! Surface pressure
!-----------------
call get_var3("ps",var3_read_1)
write(*,*) "Data for surface pressure downloaded."

! Compute average over the year for each point
var2_read = sum(var3_read_1,3)/timelen
#ifndef CPP_1D
    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,ps_timeseries)
    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,ps_avg)
#else
    ps_timeseries(1,:) = var3_read_1(1,1,:)
    ps_avg(1) = var2_read(1,1)
#endif
write(*,*) "Average of surface pressure computed."

! CO2 vmr
!--------
call get_var3("co2_layer1",var3_read_1)
where (var3_read_1 < 0.)
    var3_read_1 = 1.e-10
else where (var3_read_1 > 1.)
    var3_read_1 = 1.
end where
#ifndef CPP_1D
    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,q_co2)
#else
    q_co2(1,:) = var3_read_1(1,1,:)
#endif
A = (1./m_co2 - 1./m_noco2)
B = 1./m_noco2
vmr_co2_PCM = q_co2/(A*q_co2 + B)/m_co2
write(*,*) "Data for CO2 vmr downloaded."

! H2O vmr
!--------
call get_var3("h2o_layer1",var3_read_1)
where (var3_read_1 < 0.)
    var3_read_1 = 1.e-10
else where (var3_read_1 > 1.)
    var3_read_1 = 1.
end where
#ifndef CPP_1D
    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,q_h2o)
#else
    q_h2o(1,:) = var3_read_1(1,1,:)
#endif
write(*,*) "Data for H2O vmr downloaded."

do islope = 1,nslope
    if (nslope /= 1) then
        write(num,'(i2.2)') islope
        num = '_slope'//num
    endif

    ! CO2 ice
    !--------
    call get_var3("co2ice"//num,var3_read_1)
    where (var3_read_1 < 0.) var3_read_1 = 0.
    write(*,*) "Data for co2_ice"//num//" downloaded."
#ifndef CPP_STD
    call get_var3("perennial_co2ice"//num,var3_read_2)
    write(*,*) "Data for perennial_co2ice"//num//" downloaded."
#endif

    ! Compute the minimum over the year for each point
    var2_read = minval(var3_read_1 + var3_read_2,3)
#ifndef CPP_1D
    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,2))
#else
    min_co2_ice(1,islope,2) = var2_read(1,1)
#endif
    write(*,*) "Min of co2_ice"//num//" computed."

    ! H2O ice
    !--------
    call get_var3("h2o_ice_s"//num,var3_read_1)
    where (var3_read_1 < 0.) var3_read_1 = 0.
    write(*,*) "Data for h2o_ice_s"//num//" downloaded."
#ifndef CPP_STD
    call get_var3("watercap"//num,var3_read_2)
    write(*,*) "Data for watercap"//num//" downloaded."
#endif

    ! Compute the minimum over the year for each point
    var2_read = minval(var3_read_1 + var3_read_2,3)
#ifndef CPP_1D
    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,2))
#else
    min_h2o_ice(1,islope,2) = var2_read(1,1)
#endif
    write(*,*) "Min of h2o_ice"//num//" computed."

    ! Tsurf
    !------
    call get_var3("tsurf"//num,var3_read_1)
       write(*,*) "Data for tsurf"//num//" downloaded."

    ! Compute average over the year for each point
    var2_read = sum(var3_read_1,3)/timelen
#ifndef CPP_1D
    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg(:,islope))
#else
    tsurf_avg(1,islope) = var2_read(1,1)
#endif
    write(*,*) "Average of tsurf"//num//" computed."

    if (soil_pem) then
        ! Tsoil
        !------
        call get_var4("soiltemp"//num,var4_read)
        write(*,*) "Data for soiltemp"//num//" downloaded."

        ! Compute average over the year for each point
        var3_read_3 = sum(var4_read,4)/timelen
#ifndef CPP_1D
        do l = 1,nsoilmx
            call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),tsoil_timeseries(:,l,islope,:))
        enddo
        call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var3_read_3,tsoil_avg(:,:,islope))
#else
        tsoil_timeseries(1,:,islope,:) = var4_read(1,1,:,:)
        tsoil_avg(1,:,islope) = var3_read_3(1,1,:)
#endif
        write(*,*) "Average of tsoil"//num//" computed."

        ! Soil water density
        !-------------------
        call get_var4("waterdensity_soil"//num,var4_read)
#ifndef CPP_1D
        do l = 1,nsoilmx
            call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),watersoil_density_timeseries(:,l,islope,:))
        enddo
#else
        watersoil_density_timeseries(1,:,islope,:) = var4_read(1,1,:,:)
#endif
        write(*,*) "Data for waterdensity_soil"//num//" downloaded."

        ! Surface water density
        !----------------------
        call get_var3("waterdensity_surface"//num,var3_read_1)
        write(*,*) "Data for waterdensity_surface"//num//" downloaded."

        ! Compute average over the year for each point
        var2_read = sum(var3_read_1,3)/timelen
#ifndef CPP_1D
        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,watersurf_density_avg(:,islope))
#else
        watersurf_density_avg(1,islope) = var2_read(1,1)
#endif
        write(*,*) "Average of waterdensity_surface"//num//" computed."
    endif
enddo
deallocate(var2_read,var3_read_1,var3_read_2,var3_read_3,var4_read,num)

! Close the NetCDF file
call error_msg(nf90_close(fID),"close",filename2)
write(*,*) '"> '//filename2//'" processed!'

END SUBROUTINE read_data_PCM

!=======================================================================

SUBROUTINE check_dim(n1,n2,str1,str2)

implicit none

integer,            intent(in) :: n1, n2
character(len = *), intent(in) :: str1, str2

character(256) :: s1, s2

if (n1 /= n2) then
    s1 = 'value of '//trim(str1)//' ='
    s2 = ' read in starting file differs from parametrized '//trim(str2)//' ='
    write(msg,'(10x,a,i4,2x,a,i4)')trim(s1),n1,trim(s2),n2
    call abort_gcm(trim(modname),trim(msg),1)
endif

END SUBROUTINE check_dim

!=======================================================================

SUBROUTINE get_var1(var,v)

implicit none

character(len = *), intent(in)  :: var
real, dimension(:), intent(out) :: v

call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)

END SUBROUTINE get_var1

!=======================================================================

SUBROUTINE get_var3(var,v) ! on U grid

implicit none

character(len = *),     intent(in)  :: var
real, dimension(:,:,:), intent(out) :: v

call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)

END SUBROUTINE get_var3

!=======================================================================

SUBROUTINE get_var4(var,v)

implicit none

character(len = *),       intent(in)  :: var
real, dimension(:,:,:,:), intent(out) :: v

call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)

END SUBROUTINE get_var4

!=======================================================================

SUBROUTINE error_msg(ierr,typ,nam)

implicit none

integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME

if (ierr == nf90_noerr) return
select case(typ)
    case('inq');   msg="Field <"//trim(nam)//"> is missing"
    case('get');   msg="Reading failed for <"//trim(nam)//">"
    case('put');   msg="Writing failed for <"//trim(nam)//">"
    case('open');  msg="File opening failed for <"//trim(nam)//">"
    case('close'); msg="File closing failed for <"//trim(nam)//">"
    case default
        write(*,*) 'There is no message for this error.'
        error stop
end select
call abort_gcm(trim(modname),trim(msg),ierr)

END SUBROUTINE error_msg

END MODULE read_data_PCM_mod