PROGRAM main_ddt

    use sfx_phys_const
    USE sfx_esm_param
    USE inputdata
	USE sfx_esm

    type (meteoDataType):: data_in1
       
	type (sfxDataType) :: data_outdef1

	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

    !type :: datatype_inMAS1
    !    real, allocatable :: mas_w(:)    !
    !    real, allocatable :: mas_dt(:)
    !    real, allocatable :: mas_st(:)
    !    real, allocatable :: mas_dq(:)
    !    real, allocatable :: mas_cflh(:)
    !    real, allocatable :: mas_z0in(:)
    !end type
    type(meteoDataVecType) :: meteo

    !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

    !type :: datatype_outMAS1
    !    real, allocatable :: masout_zl(:)
    !    real, allocatable :: masout_ri(:)
    !    real, allocatable :: masout_re(:)
    !    real, allocatable :: masout_lnzuzt(:)
    !    real, allocatable :: masout_zu(:)
    !    real, allocatable :: masout_ztout(:)
    !    real, allocatable :: masout_rith(:)
    !    real, allocatable :: masout_cm(:)
    !    real, allocatable :: masout_ch(:)
    !    real, allocatable :: masout_ct(:)
    !    real, allocatable :: masout_ckt(:)
    !end type

    type(sfxDataVecType) :: data_outMAS

    !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_outMAS%zeta, data_outMAS%Rib, data_outMAS%Re, data_outMAS%B,&
            !data_outMAS%z0_m,data_outMAS%z0_t,data_outMAS%Rib_conv_lim,data_outMAS%Cm,&
            !data_outMAS%Ct,data_outMAS%Km,data_outMAS%Pr_t_inv,&
            data_par1, numst)

    !CALL surf_fluxMAS(meteo%mas_w, meteo%mas_dt, meteo%mas_st, meteo%mas_dq,&
    !        meteo%mas_cflh, meteo%mas_z0in,&
    !        data_outMAS%masout_zl, data_outMAS%masout_ri, data_outMAS%masout_re, data_outMAS%masout_lnzuzt,&
    !        data_outMAS%masout_zu,data_outMAS%masout_ztout,data_outMAS%masout_rith,data_outMAS%masout_cm,&
    !        data_outMAS%masout_ch,data_outMAS%masout_ct,data_outMAS%masout_ckt,&
    !        data_par1, data_lutyp1,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