program sfx_main

    use sfx_phys_const
    use sfx_esm_param
    use sfx_esm

    integer, parameter :: test = 1


    type(meteoDataType):: data_in1
    type(meteoDataVecType) :: meteo
	type(sfxDataType) :: data_outdef1
    type(sfxDataVecType) :: data_outMAS
    type(numericsType) :: data_par1


    integer :: numst, i
    real :: cflh, z0in
    character(len = 50) :: filename_in
    character(len = 50) :: filename_out
    character(len = 50) :: filename_in2

    !input
    !  mas_w - abs(wind velocity) at constant flux layer (cfl) hight (m/s)
    !  mas_dt - difference between potential temperature at cfl hight and at surface  ( deg. k)
    !  mas_st - semi-sum of potential temperature at cfl hight and   and at surface  ( deg. k)
    !  mas_dq - difference between humidity at cfl hight and a surface   ( gr/gr )
    !  cflh - - cfl hight ( m )
    !  z0in=0.01 - roughness of surface ( m );
    !  it    - number of iterations
    !  lu_indx -  1 for land, 2 for sea, 3 for lake
    !  test - file input


    !output
    !masout_zl - non-dimensional cfl hight
    !masout_ri - richardson number
    !masout_re  - reynods number
    !masout_lnzuzt - ln(zu/zt)
    !masout_zu  - dynamical roughness zu (m)
    !masout_ztout - thermal   roughness zt (m)
    !masout_rith - critical richardson number
    !masout_cm  - transfer coefficient for momentum
    !masout_ch - transfer coefficient fr heat
    !masout_ct - coefficient of turbulence (km) at cfl hight (m**2/s)
    !masout_ckt - alft=kt/km ( kt-coefficient of turbulence for heat)

    !> @brief Test - file selection  for test

    write(*,*) 'running code'

    if (TEST==1) then
        filename_in='data/MOSAiC.txt'
        filename_out='out_MOSAiC.txt'
        filename_in2='data/MOSAiC_zh.txt'
    elseif (TEST==2) then
        filename_in='data/Irgason1.txt'
        filename_out='out_IRGASON1.txt'
        filename_in2='data/IRGASON_zh.txt'
    endif

    open (1, file= filename_in, status ='old')
    open (2, file=filename_out)
    numst=0
    do WHILE (ioer.eq.0)
        read (1,*, iostat=ioer)  data_in1%U,  data_in1%dT,  data_in1%Tsemi,  data_in1%dQ
        numst=numst+1
    enddo
    close (1)
    numst=numst-1


    allocate(meteo%h(numst))
    allocate(meteo%U(numst))
    allocate(meteo%dT(numst))
    allocate(meteo%Tsemi(numst))
    allocate(meteo%dQ(numst))
    allocate(meteo%z0_m(numst))


    allocate(data_outMAS%zeta(numst))
    allocate(data_outMAS%Rib(numst))
    allocate(data_outMAS%Re(numst))
    allocate(data_outMAS%B(numst))
    allocate(data_outMAS%z0_m(numst))
    allocate(data_outMAS%z0_t(numst))
    allocate(data_outMAS%Rib_conv_lim(numst))
    allocate(data_outMAS%Cm(numst))
    allocate(data_outMAS%Ct(numst))
    allocate(data_outMAS%Km(numst))
    allocate(data_outMAS%Pr_t_inv(numst))

    open (11, file=filename_in2, status ='old')
    open (1, file= filename_in, status ='old')
    read (11, *) cflh, z0in
    do i=1,numst
        read (1,*) data_in1%U,  data_in1%dT,  data_in1%Tsemi,  data_in1%dQ

        meteo%h(i)=cflh
        meteo%U(i) = meteo%U(i)+data_in1%U
        meteo%dT(i) = meteo%dT(i)+data_in1%dT
        meteo%Tsemi(i) = meteo%Tsemi(i)+data_in1%Tsemi
        meteo%dQ(i) = meteo%dQ(i)+data_in1%dQ
        meteo%z0_m(i)=z0in
    enddo

    CALL get_surface_fluxes_vec(data_outMAS, meteo, &
            data_par1, numst)

    do i=1,numst
        write (2,20) data_outMAS%zeta(i), data_outMAS%Rib(i), data_outMAS%Re(i), data_outMAS%B(i),&
                data_outMAS%z0_m(i), data_outMAS%z0_t(i), data_outMAS%Rib_conv_lim(i), data_outMAS%Cm(i),&
                data_outMAS%Ct(i), data_outMAS%Km(i), data_outMAS%Pr_t_inv(i)
    enddo

    deallocate(meteo%h)
    deallocate(meteo%U)
    deallocate(meteo%dT)
    deallocate(meteo%Tsemi)
    deallocate(meteo%dQ)
    deallocate(meteo%z0_m)

    deallocate(data_outMAS%zeta)
    deallocate(data_outMAS%Rib)
    deallocate(data_outMAS%Re)
    deallocate(data_outMAS%B)
    deallocate(data_outMAS%z0_m)
    deallocate(data_outMAS%z0_t)
    deallocate(data_outMAS%Rib_conv_lim)
    deallocate(data_outMAS%Cm)
    deallocate(data_outMAS%Ct)
    deallocate(data_outMAS%Km)
    deallocate(data_outMAS%Pr_t_inv)

    10 format (f8.4,2x,f8.4)
    20 format (11(f10.4,3x))

    stop
end program