module grid

    ! интерфейс
    ! ---------------------------------------------------------------------------------

    use const,      only : miss_v, pi, r_earth, deg2rad, day2sec, len_default
    use datetime,   only : date_type
    use netcdf
    use netcdf_kit, only : nc_errhand

    implicit none

    public
    private :: spatial_grid_init
    private :: temporal_grid_init
    private :: vertical_grid_init
    
    ! переменные
    ! ---------------------------------------------------------------------------------
    
    ! широтно-долготная сетка:
    integer :: nlon, nlat
    real :: dlon, dlat, lon_ref, lat_ref
    integer :: i0           !< Longitude of the Western edge of regional domain, number of grid
    integer :: i1           !< Longitude of the Eastern edge of regional domain, number of grid
    integer :: j0           !< Latitude of the Southern edge of regional domain, number of grid
    integer :: j1           !< Latitude of the Northern edge of regional domain, number of grid
    integer :: ni, nj
    integer :: ni_nc, nj_nc, i0_nc, i1_nc, j0_nc, j1_nc
    integer :: nlon_nc, nlat_nc
    real :: dlon_nc, dlat_nc, lon_ref_nc, lat_ref_nc
    real :: lon0_nc, lon1_nc, lat0_nc, lat1_nc
    integer :: ncell_tot, ncell_mask
    integer, allocatable :: id(:,:), id_to_i(:), id_to_j(:)
    integer, allocatable :: ich(:,:), ich_to_i(:), ich_to_j(:)
    real, allocatable :: lon(:), lat(:), area(:,:)  ! SC1*SC(J)
    integer, allocatable :: mask(:,:)
    real, allocatable :: mask_nc(:,:)
    integer :: nid, nich
    real, parameter :: dlon_default = 0.5, dlat_default = 0.5

    ! вертикальная сетка:
    integer, parameter :: ms = 4                !< номер уровня по высоте, соотвествующего поверхности почвы
    integer, parameter :: ml = 27               !< номер уровня по высоте, соотвествующего нижней границе почвы
    real, parameter :: z(ml) = (/-1.,-2/3.,-1/3.,0.,1.,2.,4.,8.,15.,25.,35.,45.,55.,65.,75.,85.,95.,105.,115.,125.,135.,145.,155.,200.,300.,500.,1000./)
    real :: dz(ml)

    ! сетка по времени:
    integer :: dt            !< Timestep, sec
    integer :: ntime         !< Number of timesteps
    integer :: nt_nc
    integer :: dt_nc
    real(8) :: jday0_nc
    character(len_default) :: DL
    character(len_default) :: DI
    type(date_type) :: date_c, date_fst, date_lst
    integer :: year_min, year_max
	
    ! текущие индексы в циклах:
    integer :: tt
    integer :: ii
    integer :: jj
    integer :: nn


contains


    ! внешние процедуры
    ! ---------------------------------------------------------------------------------
    
    subroutine grid_init()
        ! ---------------------------------------
        call spatial_grid_init()
        
        call temporal_grid_init()
        
        call vertical_grid_init()
        
    end subroutine

    
    ! внутренние процедуры
    ! ---------------------------------------------------------------------------------

    subroutine spatial_grid_init()
        ! ---------------------------------------
        !< @brief широтно-долготная сетка
    
        use config, only : environment_data_type, &
                            & lsm_datafile, &
                            & nc_namespace, &
                            & spatial_grid_mode, &
                            & spatial_grid_res, &
                            & spatial_sample_mode, &
                            & id_user => id, &
                            & ich_user => ich, &
                            & point, &
                            & polygon
        use paths, only : path_inmcm_data
        
        integer :: i, j, n
        integer :: ncid, dimid, varid, ios
        real, allocatable :: work1r(:), work2r(:,:)
        integer, allocatable :: work2i(:,:)
        real :: fill_v

        if (environment_data_type == 'lsm_offline') then
                
            call nc_errhand( nf90_open(trim(lsm_datafile), nf90_nowrite, ncid) )
            
            call nc_errhand( nf90_inq_dimid(ncid, 'lon', dimid) )
            call nc_errhand( nf90_inquire_dimension(ncid, dimid, len = ni_nc) )
            call nc_errhand( nf90_inq_varid(ncid, 'lon', varid) )
            allocate(work1r(ni_nc))
            call nc_errhand( nf90_get_var(ncid, varid, work1r(:), (/1/), (/ni_nc/)) )
            if (ni_nc > 1) then
                dlon_nc = work1r(2) - work1r(1)
            else
                dlon_nc = dlon_default
            endif
            if (dlon_nc < 0.) stop 'check failed : dlon_nc < 0'
            lon0_nc = work1r(1)
            lon1_nc = work1r(ni_nc)
            deallocate(work1r)
        
            call nc_errhand( nf90_inq_dimid(ncid, 'lat', dimid) )
            call nc_errhand( nf90_inquire_dimension(ncid, dimid, len = nj_nc) )
            call nc_errhand( nf90_inq_varid(ncid, 'lat', varid) )
            allocate(work1r(nj_nc))
            call nc_errhand( nf90_get_var(ncid, varid, work1r(:), (/1/), (/nj_nc/)) )
            if (nj_nc > 1) then
                dlat_nc = work1r(2) - work1r(1)
            else
                dlat_nc = dlat_default
            endif
            if (dlat_nc < 0.) stop 'check failed : dlat_nc < 0'
            lat0_nc = work1r(1)
            lat1_nc = work1r(nj_nc)
            deallocate(work1r)
            
            nlon_nc = 360/dlon_nc
            nlat_nc = 180/dlat_nc
            lon_ref_nc = -180. + 0.5*dlon_nc
            lat_ref_nc = -90. + 0.5*dlat_nc
            
            i0_nc = nint((lon0_nc - lon_ref_nc)/dlon_nc) + 1
            i1_nc = nint((lon1_nc - lon_ref_nc)/dlon_nc) + 1
            j0_nc = nint((lat0_nc - lat_ref_nc)/dlat_nc) + 1
            j1_nc = nint((lat1_nc - lat_ref_nc)/dlat_nc) + 1
            
            allocate(mask_nc(i0_nc:i1_nc,j0_nc:j1_nc))
            ios = nf90_inq_varid(ncid, nc_namespace%mask, varid)
            if (ios == 0) then
                call nc_errhand( nf90_get_var(ncid, varid, mask_nc(:,:), (/1,1/), (/ni_nc,nj_nc/)) )
            else
                call nc_errhand( nf90_inq_varid(ncid, nc_namespace%ta, varid) )
                call nc_errhand( nf90_get_att(ncid, varid, '_FillValue', fill_v) )
                allocate(work2r(i0_nc:i1_nc,j0_nc:j1_nc))
                call nc_errhand( nf90_get_var(ncid, varid, work2r(:,:), (/1,1/), (/ni_nc,nj_nc/)) )
                do i = i0_nc, i1_nc
                    do j = j0_nc, j1_nc
                        if (work2r(i,j) == fill_v) then
                            mask_nc(i,j) = 0
                        else
                            mask_nc(i,j) = 1
                        endif
                    enddo
                enddo
            endif
            
            call nc_errhand( nf90_close(ncid) )
            
        else
        
            ni_nc = miss_v
            nj_nc = miss_v
            dlon_nc = miss_v
            dlat_nc = miss_v
            lon0_nc = miss_v
            lon1_nc = miss_v
            lat0_nc = miss_v
            lat1_nc = miss_v
            nlon_nc = miss_v
            nlat_nc = miss_v
            lon_ref_nc = miss_v
            lat_ref_nc = miss_v
            i0_nc = miss_v
            i1_nc = miss_v
            j0_nc = miss_v
            j1_nc = miss_v
            
        endif
        
        
        select case(spatial_grid_mode)
            case('none')
                dlon = miss_v
                dlat = miss_v
            case('auto')
                dlon = dlon_nc
                dlat = dlon_nc
            case('manual')
                dlon = spatial_grid_res(1)
                dlat = spatial_grid_res(2)
                if (environment_data_type == 'lsm_offline') then
                    if (dlon /= dlon_nc .or. dlat /= dlat_nc) then
                        print*, 'user grid :', dlon, dlat
                        print*, 'netcdf grid :', dlon_nc, dlat_nc
                        stop "check failed : interpolation is not supported yet"
                    endif
                endif
        end select
        
        if (spatial_grid_mode /= 'none' .and. (dlon /= 0.5 .or. dlat /= 0.5)) then
            stop "check failed : now only 05x05 grid is supported (because of INMCM external data)"
        endif
        
        select case(spatial_grid_mode)
            case('none')
                nlon = 1
                nlat = 1
                lon_ref = miss_v
                lat_ref = miss_v
            case default
                nlon = 360/dlon
                nlat = 180/dlat
                lon_ref = -180. + 0.5*dlon
                lat_ref = -90. + 0.5*dlat
        end select
                
        allocate(id(nlon,nlat))
        nid = nlon*nlat
        allocate(id_to_i(nid))
        allocate(id_to_j(nid))
        n = 0
        do i = 1, nlon
            do j = nlat, 1, -1  ! порядок обхода с СЗ на юг (как в QGIS)
                n = n + 1
                id(i,j) = n
                id_to_i(n) = i
                id_to_j(n) = j
            enddo
        enddo
        
        if (dlon == 0.5 .and. dlat == 0.5) then
            allocate(work2i(nlon+2,nlat+1))
            open(1, file=path_inmcm_data//'OLIJAN', status='old')
            read(1,'(80i1)') ((work2i(i,j),i=1,nlon+2),j=1,nlat+1)
            close(1)
            allocate(ich(nlon,nlat))
            nich = count(work2i(2:nlon+1,2:nlat+1) == 1)
            allocate(ich_to_i(nich))
            allocate(ich_to_j(nich))
            ich(:,:) = 0
            n = 0
            do j = 1, nlat
                do i = 1, nlon    ! порядок обхода с ЮЗ на восток (как в INMCM)
                    if (work2i(i+1,j+1) == 1) then
                        n = n + 1
                        ich(i,j) = n
                        ich_to_i(n) = i
                        ich_to_j(n) = j
                    endif
                enddo
            enddo
            deallocate(work2i)
        else
            allocate(ich, source = id)
            allocate(ich_to_i, source = id_to_i)
            allocate(ich_to_j, source = id_to_j)
            ich = id
            ich_to_i = id_to_i
            ich_to_j = id_to_j
            print*, 'warning : grid resolution is not suitable for OLIJAN'
            if (spatial_sample_mode == 'cell_ich') then
                stop "check failed : can not realize spatial_sample_mode = 'cell_ich'"
            endif
        endif
        
        
        select case(spatial_grid_mode)
            case('none')
            
                i0 = 1
                i1 = 1
                j0 = 1
                j1 = 1
                allocate(lon(i0:i1))
                allocate(lat(j0:j1))
                allocate(area(i0:i1,j0:j1))
                lon = miss_v
                lat = miss_v
                area = miss_v
                
            case default
            
                select case(spatial_sample_mode)
                    case('id')
                        i0 = id_to_i(id_user)
                        j0 = id_to_j(id_user)
                        i1 = i0
                        j1 = j0
                    case('ich')
                        i0 = ich_to_i(ich_user)
                        j0 = ich_to_j(ich_user)
                        i1 = i0
                        j1 = j0
                    case('point')
                        i0 = nint((point(1) - lon_ref)/dlon) + 1
                        j0 = nint((point(2) - lat_ref)/dlat) + 1
                        i1 = i0
                        j1 = j0
                    case('polygon')
                        i0 = nint((polygon(1) - lon_ref)/dlon) + 1
                        i1 = nint((polygon(2) - lon_ref)/dlon) + 1
                        j0 = nint((polygon(3) - lat_ref)/dlat) + 1
                        j1 = nint((polygon(4) - lat_ref)/dlat) + 1
                    case('all')
                        i0 = i0_nc
                        i1 = i1_nc
                        j0 = j0_nc
                        j1 = j1_nc
                end select
                
                allocate(lon(nlon))
                do i = 1, nlon
                    lon(i) = lon_ref + (i-1)*dlon
                enddo
                
                allocate(lat(nlat))
                do j = 1, nlat
                    lat(j) = lat_ref + (j-1)*dlat
                enddo
                
                allocate(area(nlon,nlat))
                do j = 1, nlat
                    do i = 1, nlon
                        area(i,j) = r_earth*abs(dlat)*deg2rad * r_earth*cos(lat(j)*deg2rad)*abs(dlon)*deg2rad
                    enddo
                enddo
                
        end select
        
        allocate(mask(i0:i1,j0:j1))
        select case(environment_data_type)
            case('lsm_offline')
                mask(:,:) = mask_nc(i0:i1,j0:j1)
            case default
                mask(:,:) = 1.
        end select
        
        ni = i1-i0+1
        nj = j1-j0+1
        ncell_tot = ni*nj
        ncell_mask = count(mask == 1)
        
    end subroutine
    


    subroutine temporal_grid_init()
        ! ---------------------------------------
        !< @brief сетка по времени
    
        use config, only : environment_data_type, &
                         & lsm_datafile, &
                         & dt_mode, &
                         & dt_user => dt, &
                         & datetime_init_mode, &
                         & timestamp_fst => datetime_init, &
                         & ntime_mode, &
                         & ntime_user => ntime, &
                         & timestamp_lst => datetime_last, &
                         & UTC_user => UTC, &
                         & datetime_init_1, &     
                         & datetime_last_1, &
                         & datetime_init_2, &
                         & datetime_last_2, &
                         & datetime_init_3, &
                         & datetime_last_3, &
                         & datetime_init_4, &
                         & datetime_last_4, &
                         & datetime_init_5, &
                         & datetime_last_5, &
                         & station_name
        
        integer :: ncid, dimid, varid
        real(8), allocatable :: work1r8(:)
        
        select case(environment_data_type)
            case('lsm_offline')
            
                call nc_errhand( nf90_open(trim(lsm_datafile), nf90_nowrite, ncid) )
                
                call nc_errhand( nf90_inq_dimid(ncid, 'time', dimid) )
                call nc_errhand( nf90_inquire_dimension(ncid, dimid, len = nt_nc) )
                call nc_errhand( nf90_inq_varid(ncid, 'time', varid) )
                allocate(work1r8(2))
                call nc_errhand( nf90_get_var(ncid, varid, work1r8(:), (/1/), (/2/)) )
                dt_nc = nint((work1r8(2) - work1r8(1)) * day2sec)
                jday0_nc = work1r8(1)
                deallocate(work1r8)
                
                call nc_errhand( nf90_close(ncid) )

            case('station')
              select case(station_name)
                case('DAO4')
                  DI = datetime_init_1
                  DL = datetime_last_1
                case('DAO3')        
                  DI = datetime_init_2
                  DL = datetime_last_2
                case('TRGK')
                  DI = datetime_init_3
                  DL = datetime_last_3
                case('VLDMR')
                  DI = datetime_init_4
                  DL = datetime_last_4
                case('Rostov')
                  DI = datetime_init_5
                  DL = datetime_last_5
              end select
               
            case default
                nt_nc = miss_v
                dt_nc = miss_v
                jday0_nc = miss_v
        end select

        select case(dt_mode)
            case('auto')
                dt = dt_nc
            case('manual')
                dt = dt_user
                if (environment_data_type == 'lsm_offline') then
                    if (dt /= dt_nc) then
                        print*, 'user dt :', dt
                        print*, 'netcdf dt :', dt_nc
                        stop "check failed : interpolation is not supported yet"
                    endif
                endif
        end select

        select case(datetime_init_mode)
            case('auto')
                select case(environment_data_type)
                  case('lsm_offline')
                     call date_fst%init(jday = jday0_nc, UTC = UTC_user)
                  case('station')
                    call date_fst%init(timestamp = DI, UTC = UTC_user)
                  case default
                    stop "check failed : this option can not work with these settings"
                end select
            case('manual')
                call date_fst%init(timestamp = timestamp_fst, UTC = UTC_user)
        end select


        select case(ntime_mode)
            case('auto')
                select case(environment_data_type)
                  case('lsm_offline')
                     ntime = nt_nc
                     call date_lst%init(jday = date_fst%jday + ntime*dble(dt)/day2sec)
                  case('station')
                    call date_lst%init(timestamp = DL, UTC = UTC_user)
                    ntime = nint((date_lst%jday - date_fst%jday) * day2sec/dt)
                  case default
                    stop "check failed : this option can not work with these settings"
                end select
            case('ntime')
                ntime = ntime_user
                call date_lst%init(jday = date_fst%jday + ntime*dble(dt)/day2sec)
            case('datetime_last')
                call date_lst%init(timestamp = timestamp_lst)
                ntime = nint((date_lst%jday - date_fst%jday) * day2sec/dt)
        end select

        date_c = date_fst

		year_min = date_fst%y
		year_max = date_lst%y

    end subroutine

    
    
    subroutine vertical_grid_init()
        ! ---------------------------------------
        !< @brief вертикальная сетка
        
        integer :: k
        
        dz(:) = miss_v
        do k = ms+1, ml-1
            dz(k) = 0.5*(z(k+1) - z(k-1))
        enddo
        
    end subroutine

        
end module