MODULE read_data_GCM_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(256) :: msg, var, modname ! for reading
integer        :: fID, vID          ! for reading

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

SUBROUTINE read_data_GCM(fichnom,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,           &
                         min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
                         watersurf_density_ave,watersoil_density)
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 GCM
!
! Authors: RV & LL
!=======================================================================

include "dimensions.h"

!=======================================================================
! Arguments:

character(len = *), intent(in) :: fichnom                             !--- FILE NAME
integer,            intent(in) :: timelen                             ! number of times stored in the file
integer                        :: 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),         intent(out) :: min_co2_ice      ! Minimum of co2 ice  per slope of the year [kg/m^2]
real, dimension(ngrid,nslope),         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_gcm_phys ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
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,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per  slope of the year [kg/m^2]
!SOIL
real, dimension(ngrid,nslope),                 intent(out) :: tsurf_ave             ! Average surface temperature of the concatenated file [K]
real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_ave             ! Average soil temperature of the concatenated file [K]
real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_gcm             ! Surface temperature of the concatenated file, time series [K]
real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_gcm             ! Soil temperature of the concatenated file, time series [K]
real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_ave ! Water density at the surface [kg/m^3]
real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
!===============================================================================
!   Local Variables
character(12)                   :: start_file_type = "earth" ! default start file type
real, dimension(:), allocatable :: time      ! times stored in start
integer                         :: indextime ! index of selected time
integer                         :: edges(4), corner(4)
integer                         :: i, j, l, t          ! loop variables
real                            ::  A , B, mmean       ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
integer                         :: islope              ! loop for variables
character(2)                    :: num                 ! for reading sloped variables
real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn         ! h2o ice per slope of the concatenated file [kg/m^2]
real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap_slope
real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm           ! CO2 volume mixing ratio in the first layer  [mol/m^3]
real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                ! Surface Pressure [Pa]
real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn         ! Average surface temperature of the concatenated file [K]
real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn         ! Average soil temperature of the concatenated file [K]
real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn         ! Surface temperature of the concatenated file, time series [K]
real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn         ! Soil temperature of the concatenated file, time series [K]
real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn             ! CO2 mass mixing ratio in the first layer [kg/m^3]
real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn             ! H2O mass mixing ratio in the first layer [kg/m^3]
real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn     ! co2 ice amount per  slope of the year [kg/m^2]
real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn ! Water density at the surface, time series [kg/m^3]
real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn ! Water density in the soil layer, time series [kg/m^3]
real, dimension(ngrid,nslope,timelen)                               :: watersurf_density     ! Water density at the surface, time series [kg/m^3]

!-----------------------------------------------------------------------
modname="read_data_gcm"

A = (1/m_co2 - 1/m_noco2)
B = 1/m_noco2

write(*,*) "Opening ", fichnom, "..."

!  Open initial state NetCDF file
var = fichnom
call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)

write(*,*) "Downloading data for vmr co2..."

call get_var3("co2_layer1"   ,q_co2_dyn)

write(*,*) "Downloading data for vmr co2 done"
write(*,*) "Downloading data for vmr h20..."

call get_var3("h2o_layer1"   ,q_h2o_dyn)

write(*,*) "Downloading data for vmr h2o done"
write(*,*) "Downloading data for surface pressure ..."

call get_var3("ps"   ,ps_GCM)

write(*,*) "Downloading data for surface pressure done"
write(*,*) "nslope=", nslope

if (nslope > 1) then
    write(*,*) "Downloading data for co2ice_slope ..."

    do islope = 1,nslope
        write(num,'(i2.2)') islope
        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
    enddo
    
    write(*,*) "Downloading data for co2ice_slope done"
    write(*,*) "Downloading data for h2o_ice_s_slope ..."
    
    do islope = 1,nslope
        write(num,'(i2.2)') islope
        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
    enddo

    write(*,*) "Downloading data for h2o_ice_s_slope done"

#ifndef CPP_STD
    write(*,*) "Downloading data for watercap_slope ..."
    do islope = 1,nslope
        write(num,'(i2.2)') islope
        call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
!        watercap_slope(:,:,:,:) = 0.
    enddo
    write(*,*) "Downloading data for watercap_slope done"
#endif

    write(*,*) "Downloading data for tsurf_slope ..."

    do islope=1,nslope
        write(num,'(i2.2)') islope
        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
    enddo

    write(*,*) "Downloading data for tsurf_slope done"

#ifndef CPP_STD
    if (soil_pem) then
        write(*,*) "Downloading data for soiltemp_slope ..."

        do islope = 1,nslope
            write(num,'(i2.2)') islope
            call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
        enddo

        write(*,*) "Downloading data for soiltemp_slope done"
        write(*,*) "Downloading data for watersoil_density ..."

        do islope = 1,nslope
            write(num,'(i2.2)') islope
            call get_var4("Waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
        enddo

        write(*,*) "Downloading data for  watersoil_density  done"
        write(*,*) "Downloading data for  watersurf_density  ..."

        do islope = 1,nslope
            write(num,'(i2.2)') islope
            call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
        enddo

        write(*,*) "Downloading data for  watersurf_density  done"

    endif !soil_pem
#endif
else !nslope = 1 no slope, we copy all the values
    call get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
    call get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
    call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
    
#ifndef CPP_STD
!    call get_var3("watercap", watercap_slope(:,:,1,:))
    watercap_slope(:,:,1,:) = 0.
    watersurf_density_dyn(:,:,:,:) = 0.
    watersoil_density_dyn(:,:,:,:,:) = 0.
#endif

    if (soil_pem) call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
endif !nslope=1

! Compute the minimum over the year for each point
write(*,*) "Computing the min of h2o_ice_slope"
min_h2o_ice_dyn(:,:,:) = minval(h2o_ice_s_dyn + watercap_slope,4)
write(*,*) "Computing the min of co2_ice_slope"
min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn,4)

!Compute averages
write(*,*) "Computing average of tsurf"
tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen

#ifndef CPP_1D
    do islope = 1,nslope
        do t = 1,timelen
            call gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
        enddo
    enddo
#endif

if (soil_pem) then
    write(*,*) "Computing average of tsoil"
    tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
    write(*,*) "Computing average of watersurf_density"
    watersurf_density_ave(:,:) = sum(watersurf_density(:,:,:),3)/timelen
endif

! By definition, a density is positive, we get rid of the negative values
where (min_co2_ice_dyn < 0.) min_co2_ice_dyn = 0.
where (min_h2o_ice_dyn < 0.) min_h2o_ice_dyn = 0.

do i = 1,iim + 1
    do j = 1,jjm_input + 1
        do t = 1, timelen
            if (q_co2_dyn(i,j,t) < 0) then
                q_co2_dyn(i,j,t) = 1.e-10
            else if (q_co2_dyn(i,j,t) > 1) then
                q_co2_dyn(i,j,t) = 1.
            endif
            if (q_h2o_dyn(i,j,t) < 0) then
                q_h2o_dyn(i,j,t) = 1.e-30
            else if (q_h2o_dyn(i,j,t) > 1) then
                q_h2o_dyn(i,j,t) = 1.
            endif
            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
            vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
        enddo
    enddo
enddo

#ifndef CPP_1D
    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_GCM,ps_timeseries)
    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)

    do islope = 1,nslope
        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
        if (soil_pem) then
            do l = 1,nsoilmx
                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
                do t=1,timelen
                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
                enddo
            enddo
        endif !soil_pem
        do t = 1,timelen
            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
        enddo
    enddo

    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_ave_dyn,tsurf_ave)
#else
    vmr_co2_gcm_phys(1,:) = vmr_co2_gcm(1,1,:)
    ps_timeseries(1,:) = ps_GCM(1,1,:)
    q_co2(1,:) = q_co2_dyn(1,1,:)
    q_h2o(1,:) = q_h2o_dyn(1,1,:)
    min_co2_ice(1,:) = min_co2_ice_dyn(1,1,:)
    min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
    if (soil_pem) then
        tsoil_ave(1,:,:) = tsoil_ave_dyn(1,1,:,:)
        tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
    endif !soil_pem
    tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)
    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
    tsurf_ave(1,:) = tsurf_ave_dyn(1,1,:)
#endif

END SUBROUTINE read_data_gcm

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('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_GCM_mod