! Created by Andrey Debolskiy on 12.10.2024.
program gabls1
    use ISO_C_BINDING, only: C_NULL_CHAR
    use parkinds, only: rf=>kind_rf, im=>kind_im
    use phys_fluid
    use scm_io_default
#ifdef USE_NETCDF
    use io, only: variable_metadata, write_2d_real_nc, write_1d_real_nc
#endif
    use scm_sfx_data, only: taux, tauy, umst, hflx, qflx, cu
    use scm_state_data
    use pbl_turb_data
    use pbl_turb_common
    use pbl_solver
    use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z, &
            get_sigma_from_z_theta, get_theta_v, get_density_theta_vector
    use diag_pbl
    use pbl_grid
    use sfx_data, only: meteoDataType, sfxDataType
    use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
            numericsType_most => numericsType
    use config_parser
    use config_utils


    implicit none
    !local varaibles
    real time_begin, time_end, time_current
    real dt, out_dt, output_dt
    integer status, ierr
    logical bstatus
    character(len = 160) ::fname
    character(len = 128) :: fname_nc
    character(len = 128) ::fname_config = '../config-examples/config-gabls1-ex.txt'
    real, allocatable :: ttemp(:,:), utemp(:,:),vtemp(:,:), ttime(:)
    real(kind=rf):: ug, vg
    real(kind=rf):: z0, zh, f_cor
    real(kind=rf):: tsurf
    real(kind=rf):: shir
    real(kind=rf):: hbl_0
    real(kind=rf):: tgrad
    real(kind=rf):: cooling_rate
    type(fluidParamsDataType) fluid_params
    type(pblgridDataType) grid

    type(stateBLDataType):: bl, bl_old;
    type(turbBLDataType):: turb

    type(meteoDataType) :: meteo_cell
    type(sfxDataType) :: sfx_cell
    type(numericsType_most) :: numerics_sfx


    !diagnostic variables
    real hpbl
    integer nkpbl

    !io variables
    real, dimension (5):: series_val
    type (io_struct) :: series_f, scan_f
#ifdef USE_NETCDF
    type(variable_metadata) :: meta_z, meta_t
    type(variable_metadata) :: meta_temperature, meta_uwind, meta_vwind
#endif
    integer k, nt, kmax
    ! init experiment params
    numerics_sfx%maxiters_charnock = 10

    ug = 8.00000
    vg = 0.0
    dt = 1.0;
!    call set_fluid_default(fluid_params)
!    fluid_params%tref = 265.0
!    fluid_params%pref = 1013.2500
!    fluid_params%beta = fluid_params%g / fluid_params%tref
!    fluid_params%kappa = 0.4
!    fluid_params%p0= fluid_params%pref



    !set output filenames
    fname = 'test.dat'
    fname_nc = 'gabls1.nc'
    series_f%fname = 'phys.dat'
    call set_file(series_f, series_f%fname)

    scan_f%fname='time_scan.dat'
    call set_file(scan_f, scan_f%fname)
#ifdef USE_NETCDF
    call set_meta_nc(meta_z, meta_t, meta_temperature, meta_uwind, meta_vwind)
#endif



    call init_config(fname_config,status, ierr)
    if (status == 0) then
        write(*, *) ' FAILURE! > could not initialize config file: ', fname_config
        ierr = 1        ! signal ERROR
        stop
    end if

    call c_config_is_varname("time.begin"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_float("time.begin"//C_NULL_CHAR, time_begin, status)
        if (status == 0) then
            write(*, *) ' FAILURE! > could not set:  time.begin '
            ierr = 1        ! signal ERROR
            stop
        end if
    end if
    call c_config_is_varname("time.end"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_float("time.end"//C_NULL_CHAR, time_end, status)
        if (status == 0) then
            write(*, *) ' FAILURE! > could not set:  time.end '
            ierr = 1        ! signal ERROR
            stop
        end if
    end if
    call c_config_is_varname("time.dt"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_float("time.dt"//C_NULL_CHAR, dt, status)
        if (status == 0) then
            write(*, *) ' FAILURE! > could not set:  time.dt '
            ierr = 1        ! signal ERROR
            stop
        end if
    end if

    call c_config_is_varname("time.out_dt"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_float("output.dt"//C_NULL_CHAR, output_dt, status)
        if (status == 0) then
            write(*, *) ' FAILURE! > could not set:  output.dt '
            ierr = 1        ! signal ERROR
            stop
        end if
    end if

    call c_config_is_varname("Ug"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_float("Ug"//C_NULL_CHAR, ug, status)
        if (status == 0) then
            write(*, *) ' FAILURE! > could not set:  Ug '
            ierr = 1        ! signal ERROR
            stop
        end if
    end if
    call c_config_is_varname("Vg"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_float("Vg"//C_NULL_CHAR, vg, status)
        if (status == 0) then
            write(*, *) ' FAILURE! > could not set:  Vg '
            ierr = 1        ! signal ERROR
            stop
        end if
    end if
    call c_config_is_varname("hbl_0"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_float("hbl_0"//C_NULL_CHAR, hbl_0, status)
        if (status == 0) then
            write(*, *) ' FAILURE! > could not set:  hbl_0 '
            ierr = 1        ! signal ERROR
            stop
        end if
    end if
    call c_config_is_varname("Tgrad"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_float("Tgrad"//C_NULL_CHAR, tgrad, status)
        if (status == 0) then
            write(*, *) ' FAILURE! > could not set:  Tgrad '
            ierr = 1        ! signal ERROR
            stop
        end if
    end if
    call set_fluid_default(fluid_params)
    call get_fluid_params(fluid_params, status, ierr)
    if (status == 0) then
        write(*, *) ' FAILURE! > could not set:  Fluid params '
        ierr = 1        ! signal ERROR
        stop
    end if
    call c_config_is_varname("cooling_rate"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_float("cooling_rate"//C_NULL_CHAR, cooling_rate, status)
        if (status == 0) then
            write(*, *) ' FAILURE! > could not set:  cooling_rate '
            ierr = 1        ! signal ERROR
            stop
        end if
    end if
    call c_config_is_varname("fcoriolis"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_float("fcoriolis"//C_NULL_CHAR, f_cor, status)
        if (status == 0) then
            write(*, *) ' FAILURE! > could not set:  fcoriolis '
            ierr = 1        ! signal ERROR
            stop
        end if
    end if
    call set_fluid_default(fluid_params)
    call get_fluid_params(fluid_params, status, ierr)
    if (status == 0) then
        write(*, *) ' FAILURE! > could not set:  Fluid params '
        ierr = 1        ! signal ERROR
        stop
    end if
    !z coordinate
    call get_grid_params(grid, status, ierr)
    if (status == 0) then
        write(*, *) ' FAILURE! > could not set:  grid '
        ierr = 1        ! signal ERROR
        stop
    end if

    kmax = grid%kmax
    time_current = time_begin
    ! hack to insure shir array is good
    !shir = 3.141592653589793238 * 72.3798 / 180.0
    !f_cor = 2.0 * fluid_params%omega * sin(shir)
    bl_old%cor= f_cor
    bl%cor = f_cor
    bl%land_type = 0.0
    bl%p0=fluid_params%p0


    !setup surface
    z0 = 0.01
    zh = z0 / 10.0
    ! init blData states

    call scm_data_allocate(bl, kmax)
    call scm_data_allocate(bl_old, kmax)
    call init_state(bl, ug, vg, fluid_params%tref)
    call scm_data_copy(bl_old,bl)




    DO K = 1, kmax
        if (grid%z_cell(k) > 100.0) then
            bl_old%theta(k) = fluid_params%tref + tgrad * (grid%z_cell(k) - 100.0)
        else
            bl_old%theta(k) = fluid_params%tref
        end if
        !write(*,*) k, sig(k)* fluid_params%pref* 100.0, grid%z_cell(k)
    END DO
    !finish updating grid
    call get_sigma_from_z_theta( grid, bl, fluid_params%p0 * 100.0, fluid_params)
    call theta2ta(bl_old%temp, bl_old%theta, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
    call theta2ta(bl%temp, bl%theta, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
    call get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
    call scm_data_copy(bl_old,bl)

    !Initialize turbulent closure
    call turb_data_allocate(turb, kmax)
    ! getclosurefromconfig
    call get_closure_params(turb, status, ierr)




    !io setup
    nt = floor((time_end - time_begin)/output_dt)
    allocate(ttime(nt))
    allocate(ttemp(kmax,nt),utemp(kmax,nt),vtemp(kmax,nt))
#ifdef USE_NETCDF
    call write_1d_real_nc(grid%z_cell,fname_nc,meta_z)
#endif
    out_dt = 0.0
    nt=0
    write(*,*)'nt=0', nt
    do while (time_current < time_end)

        tsurf = fluid_params%tref - cooling_rate * (time_current - time_begin) / 3600.0;
        meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%u(kmax)**2)
        meteo_cell%dT = -tsurf + bl_old%theta(kmax)
        meteo_cell%Tsemi = 0.5 * (tsurf + bl_old%theta(kmax))
        meteo_cell%dQ  = 0.0
        meteo_cell%h = grid%z_cell(kmax)
        meteo_cell%z0_m = z0


        call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
        bl%surf%cm2u = sfx_cell%Cm**2 * meteo_cell%U
        taux = 1.1 * sfx_cell%Cm**2 * meteo_cell%U *  bl_old%u(kmax)
        tauy = 1.1 * sfx_cell%Cm**2 * meteo_cell%U *  bl_old%v(kmax)
        hflx =  sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U * meteo_cell%dT
        bl%surf%hs = bl%rho(grid%kmax) * hflx
        bl%surf%es = sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U *meteo_cell%dQ
        turb%ustr = sfx_cell%Cm * meteo_cell%U
        umst = turb%ustr

        call  get_density_theta_vector( bl, fluid_params, grid)
        call  get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
        call  get_coeffs_general(turb, bl, fluid_params, 1,  grid)
        write(*,*) 'bl%kpbl', bl%kpbl
        call solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
        call add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
        call scm_data_copy(bl_old,bl)


        time_current = time_current + dt
        if (time_current >=out_dt) then
            nt = nt+1
            ttime(nt) = time_current/3600.0
            out_dt = out_dt + output_dt
            write(*,*)'nt= ', nt


            call ta2theta( bl_old%theta, bl_old%temp, &
                    fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)

            call diag_pblh_inmcm(grid%z_cell,bl_old%theta,shir,0.0,umst,bl_old%kmax,hpbl)
            !create output
            series_val(1) = time_current/3600.0
            series_val(2) = hflx
            series_val(3) = turb%ustr
            series_val(4) = hpbl
            series_val(5) = 0.0

            call write_series(series_val, 5, series_f)

            call put_tscan( time_current/3600.0, grid%z_cell,bl_old, bl_old%kmax, scan_f)
#ifdef USE_NETCDF
            do k = 1, kl
                ttemp(k,nt) = thold(k)
                utemp(k,nt) = uold(k)
                vtemp(k,nt) = vold(k)
            end do
#endif
            !write(*,*) ccld * 10000.0
        end if
    end do
#ifdef USE_NETCDF
    write(*,*)'nt= ', nt
    call write_1d_real_nc(ttime,fname_nc,meta_t)
    write(*,*)'here'
    call  write_2d_real_nc(ttemp,fname_nc,meta_temperature)
    call  write_2d_real_nc(utemp,fname_nc,meta_uwind)
    call  write_2d_real_nc(vtemp,fname_nc,meta_vwind)
#endif

    call ta2theta(bl_old%theta, bl_old%temp, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
    call to_file_1d_3var(fname, grid%z_cell, bl_old%u, bl_old%v, bl_old%kmax)
    call to_file_1d_3var('coeffs.dat',grid%z_cell, bl_old%km, bl_old%kh, bl_old%kmax)
    call to_file_1d_2var('temperature.dat', grid%z_cell, bl_old%theta, bl_old%kmax)
    call close_file(series_f)
    call close_file(scan_f)

    call deallocate_pbl_grid(grid)
    call scm_data_deallocate(bl)
    call scm_data_deallocate(bl_old)
    call pbl_data_deallocate(turb)
end program gabls1


subroutine init_state(bl, ug_,vg_,tref)
    use parkinds, only: rf=>kind_rf, im=>kind_im
    use scm_state_data, only:stateBLDataType

    implicit none
    real(kind=rf):: ug_,vg_,tref
    type(stateBLDataType), intent(inout) :: bl
    integer kmax
    kmax= bl%kmax
    write(*,*) 'blkmax=', kmax
    bl%u(1:bl%kmax) = ug_
    bl%v(1:bl%kmax) = vg_
    bl%temp(1:bl%kmax) = tref
    bl%theta(1:bl%kmax) = tref
    bl%qv(1:bl%kmax) = 0.0
    bl%lwp(1:bl%kmax) = 0.0
    bl%cloud_frac(1:bl%kmax) = 0.0
    bl%s_e(1:bl%kmax) = 0.0
    bl%km(1:bl%kmax) = 0.0
    bl%kh(1:bl%kmax) = 0.0
    bl%vdcuv(1:bl%kmax) = 0.0
    bl%vdctq(1:bl%kmax) = 0.0
end subroutine init_state



subroutine add_coriolis(bl, bl_old, ug, vg, f, dt, grid)
    use scm_state_data
    use pbl_grid, only: pblgridDataType
    implicit none
    real, intent(in) :: ug, vg, f, dt
    type(stateBLDataType), intent(inout):: bl
    type(stateBLDataType), intent(in):: bl_old
    type(pblgridDataType), intent(in):: grid

    integer k
        do k = 1, grid%kmax
        bl%v(k) = bl_old%v(k)  - dt * f * (bl_old%u(k) - ug)
        bl%u(k) = bl_old%u(k) + dt * f * ( bl_old%v(k) - vg)
    end do

end subroutine add_coriolis

subroutine put_tscan( time, z, bl, nl, f)
    use parkinds, only : rf=>kind_rf, im => kind_im
    use scm_io_default
    use scm_state_data, only:stateBLDataType
    implicit none
    type(stateBLDataType), intent(in):: bl
    real(kind=rf), intent(in), dimension(nl):: z
    real(kind=rf),intent(in) :: time
    type (io_struct),intent(in) :: f
    integer(kind=im), intent(in) :: nl


    !local
    real(kind=rf), dimension(5,nl)::stamp
    integer k

    ! copy to stamp
    do k=1,nl
        stamp(1,k) = time
        stamp(2,k) = z(k)
        stamp(3,k) = bl%theta(k)
        stamp(4,k) = bl%u(k)
        stamp(5,k) = bl%v(k)
    end do

    ! call to write timestamp
    call write_timescan(stamp,nl, 5, f)

end subroutine put_tscan

subroutine get_surface_from_config(model, surface, z0,zh )
    use iso_c_binding, only: C_NULL_CHAR
    use config_parser
    use sfx_common, only: char_array2str
    use sfx_config
    use sfx_surface, only: get_surface_id
    use parkinds, only: rf=>kind_rf
    implicit none
    character, allocatable :: config_field(:)

    integer,intent(out) :: model, surface
    real(kind=rf), intent(out) :: z0,zh

    integer status, ierr


    call c_config_is_varname("model.id"//C_NULL_CHAR, status)
    if (status /= 0) then
        call c_config_get_string("model.id"//C_NULL_CHAR, config_field, status)
        if (status == 0) then
            ierr = 1        ! signal ERROR
            stop
        end if
        model = get_model_id(char_array2str(config_field))
        if (model == -1) then
            write(*, *) ' FAILURE! > unknown model [key]: ', trim(char_array2str(config_field))
            ierr = 1        ! signal ERROR
            return
        end if
    end if

    call c_config_is_varname("surface.type"//C_NULL_CHAR, status)
    if ((status /= 0)) then
        !< mandatory in user dataset
        call c_config_get_string("surface.type"//C_NULL_CHAR, config_field, status)
        if (status == 0) then
            ierr = 1        ! signal ERROR
            return
        end if
        surface = get_surface_id(char_array2str(config_field))
        if (surface == -1) then
            write(*, *) ' FAILURE! > unknown surface [key]: ', trim(char_array2str(config_field))
            ierr = 1        ! signal ERROR
            return
        end if
    endif


end subroutine get_surface_from_config
#ifdef USE_NETCDF
subroutine set_meta_nc(z_meta, t_meta, theta_meta, uwind_meta, vwind_meta)
    use io, only: variable_metadata
    type(variable_metadata), intent(inout):: z_meta, t_meta, theta_meta, uwind_meta, vwind_meta
    z_meta = variable_metadata( &
            name = 'Z', &
            dim_names = [character(len=32) :: 'Z', '', '', ''], &
            char_name = [character(len=32) :: &
                    'long_name', &
                    'units', &
                    '', '', '', '', '', '', '', ''], &
            char_value = [character(len=32) :: &
                    'Height', &
                    'm', &
                    '', '', '', '', '', '', '', ''], &
            real_name = [character(len=32) :: &
                    'missing_value', &
                    'scale_factor', &
                    '', '', '', '', '', '', '', ''], &
            real_value = [&
                    -9999.9, &
                            1.0, &
                            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    t_meta = variable_metadata( &
            name = 'Time', &
            dim_names = [character(len=32) :: 'Time', '', '', ''], &
            char_name = [character(len=32) :: &
                    'long_name', &
                    'units', &
                    '', '', '', '', '', '', '', ''], &
            char_value = [character(len=32) :: &
                    'Time', &
                    'hours since 00-00-00', &
                    '', '', '', '', '', '', '', ''], &
            real_name = [character(len=32) :: &
                    'missing_value', &
                    'scale_factor', &
                    '', '', '', '', '', '', '', ''], &
            real_value = [&
                    -9999.9, &
                            1.0, &
                            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    theta_meta = variable_metadata( &
            name = 'th', &
            dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
            char_name = [character(len=32) :: &
                    'long_name', &
                    'units', &
                    '', '', '', '', '', '', '', ''], &
            char_value = [character(len=32) :: &
                    'Potential Temperature', &
                    'K', &
                    '', '', '', '', '', '', '', ''], &
            real_name = [character(len=32) :: &
                    'missing_value', &
                    'scale_factor', &
                    '', '', '', '', '', '', '', ''], &
            real_value = [&
                    -9999.9, &
                            1.0, &
                            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    uwind_meta = variable_metadata( &
            name = 'ua', &
            dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
            char_name = [character(len=32) :: &
                    'long_name', &
                    'units', &
                    '', '', '', '', '', '', '', ''], &
            char_value = [character(len=32) :: &
                    'Longitude wind component', &
                    'm/s', &
                    '', '', '', '', '', '', '', ''], &
            real_name = [character(len=32) :: &
                    'missing_value', &
                    'scale_factor', &
                    '', '', '', '', '', '', '', ''], &
            real_value = [&
                    -9999.9, &
                            1.0, &
                            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    vwind_meta = variable_metadata( &
            name = 'va', &
            dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
            char_name = [character(len=32) :: &
                    'long_name', &
                    'units', &
                    '', '', '', '', '', '', '', ''], &
            char_value = [character(len=32) :: &
                    'Latitude wind component', &
                    'm/s', &
                    '', '', '', '', '', '', '', ''], &
            real_name = [character(len=32) :: &
                    'missing_value', &
                    'scale_factor', &
                    '', '', '', '', '', '', '', ''], &
            real_value = [&
                    -9999.9, &
                            1.0, &
                            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
end subroutine set_meta_nc
#endif