    MODULE drag3
        USE param

        implicit none

        type, public :: meteoDataType
            real :: h       ! constant flux layer height [m]
            real :: U       ! abs(wind speed) at 'h' [m/s]
            real :: dT      ! difference between potential temperature at 'h' and at surface [K]
            real :: Tsemi   ! semi-sum of potential temperature at 'h' and at surface [K]
            real :: dQ      ! difference between humidity at 'h' and at surface [g/g]
            real :: z0_m    ! surface aerodynamic roughness (should be < 0 for water bodies surface)
        end type

        type, public :: meteoDataVecType
            real, allocatable :: h(:)       ! constant flux layer height [m]
            real, allocatable :: U(:)       ! abs(wind speed) at 'h' [m/s]
            real, allocatable :: dT(:)      ! difference between potential temperature at 'h' and at surface [K]
            real, allocatable :: Tsemi(:)   ! semi-sum of potential temperature at 'h' and at surface [K]
            real, allocatable :: dQ(:)      ! difference between humidity at 'h' and at surface [g/g]
            real, allocatable :: z0_m(:)    ! surface aerodynamic roughness (should be < 0 for water bodies surface)
        end type

        type, public :: sfxDataType
            real :: zeta            ! = z/L [n/d]
            real :: Rib             ! bulk Richardson number [n/d]
            real :: Re              ! Reynolds number [n/d]
            real :: lnzuzt
            real :: z0_m            ! aerodynamic roughness [m]
            real :: z0_t            ! thermal roughness [m]
            real :: Rib_conv_lim    ! Ri-bulk convection critical value [n/d]
            real :: cm
            real :: ch
            real :: ct
            real :: ckt
        end type

        type, public :: numericsType
            integer :: maxiters_convection = 10    ! maximum (actual) number of iterations in convection
            integer :: maxiters_charnock = 10      ! maximum (actual) number of iterations in charnock roughness
        end type

        type, public:: data_lutyp
            integer, public :: lu_indx=1
        end type

    contains

        subroutine surf_flux_vec(meteo, &
                masout_zl,  masout_ri, masout_re, masout_lnzuzt, masout_zu, masout_ztout,&
                masout_rith, masout_cm, masout_ch, masout_ct, masout_ckt,&
                numerics, lu1,numst)

            type (meteoDataVecType), intent(in) :: meteo
            type (numericsType), intent(in) :: numerics

            real, dimension (numst) :: masout_zl
            real, dimension (numst) :: masout_ri
            real, dimension (numst) :: masout_re
            real, dimension (numst) :: masout_lnzuzt
            real, dimension (numst) :: masout_zu
            real, dimension (numst) :: masout_ztout
            real, dimension (numst) :: masout_rith
            real, dimension (numst) :: masout_cm
            real, dimension (numst) :: masout_ch
            real, dimension (numst) :: masout_ct
            real, dimension (numst) :: masout_ckt

            type (meteoDataType)  meteo_cell
            type (sfxDataType) sfx_cell

            type (data_lutyp) lu1
            integer i,numst

            do i = 1, numst
                meteo_cell = meteoDataType(&
                        h = meteo%h(i), &
                        U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
                        z0_m = meteo%z0_m(i))

                call surf_flux(meteo_cell, sfx_cell, numerics, lu1)

                masout_zl(i)=sfx_cell%zeta
                masout_ri(i)=sfx_cell%Rib
                masout_re(i)=sfx_cell%Re
                masout_lnzuzt(i)=sfx_cell%lnzuzt
                masout_zu(i)=sfx_cell%z0_m
                masout_ztout(i)=sfx_cell%z0_t
                masout_rith(i)=sfx_cell%Rib_conv_lim
                masout_cm(i)=sfx_cell%cm
                masout_ch(i)=sfx_cell%ch
                masout_ct(i)=sfx_cell%ct
                masout_ckt(i)=sfx_cell%ckt
            enddo
        END SUBROUTINE  surf_flux_vec

        function f_m_conv(zeta)
            real f_m_conv
            real zeta

            f_m_conv = (1.0 - alpha_m * zeta)**0.25

        end function f_m_conv

        function f_h_conv(zeta)
            real f_h_conv
            real zeta

            f_h_conv = (1.0 - alpha_h * zeta)**0.5

        end function f_h_conv


        ! stability functions (neutral)
        ! --------------------------------------------------------------------------------
        subroutine get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
            ! ----------------------------------------------------------------------------
            real, intent(out) :: psi_m, psi_h, zeta
            real, intent(in) :: h0_m, h0_t
            real, intent(in) :: B
            ! ----------------------------------------------------------------------------

            zeta = 0.0
            psi_m = log(h0_m)
            psi_h = log(h0_t) / Pr_t_0_inv
            if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv

        end subroutine
        ! --------------------------------------------------------------------------------

        ! stability functions (stable)
        ! --------------------------------------------------------------------------------
        subroutine get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
            ! ----------------------------------------------------------------------------
            real, intent(out) :: psi_m, psi_h, zeta
            real, intent(in) :: Rib
            real, intent(in) :: h0_m, h0_t
            real, intent(in) :: B

            real :: Rib_coeff
            real :: psi0_m, psi0_h
            real :: phi
            real :: c
            ! ----------------------------------------------------------------------------

            psi0_m = alog(h0_m)
            psi0_h = B / psi0_m

            Rib_coeff = beta_m * Rib
            c = (psi0_h + 1.0) / Pr_t_0_inv - 2.0 * Rib_coeff
            zeta = psi0_m * (sqrt(c**2 + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / &
                    (2.0 * beta_m * (1.0 - Rib_coeff))

            phi = beta_m * zeta

            psi_m = psi0_m + phi
            psi_h = (psi0_m + B) / Pr_t_0_inv + phi

        end subroutine
        ! --------------------------------------------------------------------------------

        ! stability functions (convection-semistrong)
        ! --------------------------------------------------------------------------------
        subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)
            ! ----------------------------------------------------------------------------
            real, intent(out) :: psi_m, psi_h, zeta
            real, intent(in) :: Rib
            real, intent(in) :: h0_m, h0_t
            real, intent(in) :: B
            integer, intent(in) :: maxiters

            real :: zeta0_m, zeta0_h
            real :: x0, x1, y0, y1

            integer :: i
            ! ----------------------------------------------------------------------------

            psi_m = log(h0_m)
            psi_h = log(h0_t)
            if (abs(B) < 1.0e-10) psi_h = psi_m

            zeta = Rib * Pr_t_0_inv * psi_m**2 / psi_h

            do i = 1, maxiters + 1
                zeta0_m = zeta / h0_m
                zeta0_h = zeta / h0_t
                if (abs(B) < 1.0e-10) zeta0_h = zeta0_m

                y1 = (1.0 - alpha_m * zeta)**0.25e0
                x1 = sqrt(1.0 - alpha_h_fix * zeta)

                y0 = (1.0 - alpha_m * zeta0_m)**0.25e0
                x0 = sqrt(1.0 - alpha_h_fix * zeta0_h)

                y0 = max(y0, 1.000001e0)
                x0 = max(x0, 1.000001e0)

                psi_m = log((y1 - 1.0e0)*(y0 + 1.0e0)/((y1 + 1.0e0)*(y0 - 1.0e0))) + 2.0e0*(atan(y1) - atan(y0))
                psi_h = log((x1 - 1.0e0)*(x0 + 1.0e0)/((x1 + 1.0e0)*(x0 - 1.0e0))) / Pr_t_0_inv

                if (i == maxiters + 1) exit

                zeta = Rib * psi_m**2 / psi_h
            end do

        end subroutine
        ! --------------------------------------------------------------------------------

        ! stability functions (convection)
        ! --------------------------------------------------------------------------------
        subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, &
                zeta_conv_lim, x10, y10, &
                h0_m, h0_t, B, maxiters)
            ! ----------------------------------------------------------------------------
            real, intent(out) :: psi_m, psi_h, zeta
            real, intent(in) :: Rib
            real, intent(in) :: zeta_conv_lim
            real, intent(in) :: x10, y10
            real, intent(in) :: h0_m, h0_t
            real, intent(in) :: B
            integer, intent(in) :: maxiters

            real :: zeta0_m, zeta0_h
            real :: x0, y0
            real :: p1, p0
            real :: a1, b1, c, f

            integer :: i
            ! ----------------------------------------------------------------------------

            p1 = 2.0 * atan(y10) + log((y10 - 1.0) / (y10 + 1.0))
            p0 = log((x10 - 1.0) / (x10 + 1.0))

            zeta = zeta_conv_lim
            do i = 1, maxiters + 1
                zeta0_m = zeta / h0_m
                zeta0_h = zeta / h0_t
                if (abs(B) < 1.0e-10) zeta0_h = zeta0_m

                a1 = (zeta_conv_lim / zeta)**(1.0 / 3.0)
                x0 = sqrt(1.0 - alpha_h_fix * zeta0_h)
                y0 = (1.0 - alpha_m * zeta0_m)**0.25
                c = log((x0 + 1.0)/(x0 - 1.0))
                b1 = -2.0*atan(y0) + log((y0 + 1.0)/(y0 - 1.0))
                f = 3.0 * (1.0 - a1)

                psi_m = f / y10 + p1 + b1
                psi_h = (f / x10 + p0 + c) / Pr_t_0_inv

                if (i == maxiters + 1) exit
                zeta = Rib * psi_m**2 / psi_h
            end do

        end subroutine
        ! --------------------------------------------------------------------------------

        ! define convection limit
        ! --------------------------------------------------------------------------------
        subroutine get_convection_lim(zeta_lim, Rib_lim, x10, y10, &
                h0_m, h0_t, B)
            ! ----------------------------------------------------------------------------
            real, intent(out) :: zeta_lim, Rib_lim
            real, intent(out) :: x10, y10

            real, intent(in) :: h0_m, h0_t
            real, intent(in) :: B

            real :: an4, an5
            real :: psi_m, psi_h
            ! ----------------------------------------------------------------------------

            an5=(Pr_t_inf_inv / Pr_t_0_inv)**4
            zeta_lim = (2.0E0*alpha_h-an5*alpha_m-SQRT((an5*alpha_m)**2+4.0E0*an5*alpha_h*(alpha_h-alpha_m))) &
                    / (2.0E0*alpha_h**2)

            y10 = f_m_conv(zeta_lim)
            x10 = f_h_conv(zeta_lim)

            !     ......definition of r-prim......
            an4 = zeta_lim / h0_m
            an5 = zeta_lim / h0_t
            if (abs(B) < 1.0e-10) an5=an4

            an5 = sqrt(1.0 - alpha_h_fix * an5)
            an4 = (1.0 - alpha_m * an4)**0.25

            psi_m = 2.0 * (atan(y10) - atan(an4)) + alog((y10 - 1.0) * (an4 + 1.0)/((y10 + 1.0) * (an4 - 1.0)))
            psi_h = alog((x10 - 1.0) * (an5 + 1.0)/((x10 + 1.0) * (an5 - 1.0))) / Pr_t_0_inv

            Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)

        end subroutine
        ! --------------------------------------------------------------------------------

        ! charnock roughness
        ! --------------------------------------------------------------------------------
        subroutine get_charnock_roughness(z0_m, u_dyn0, U, h, maxiters)
            ! ----------------------------------------------------------------------------
            real, intent(out) :: z0_m, u_dyn0

            real, intent(in) :: U, h
            integer, intent(in) :: maxiters

            real :: U10m
            real :: a1, c1, y1, c1min, f

            integer :: i, j
            ! ----------------------------------------------------------------------------

            U10m = U
            a1 = 0.0e0
            y1 = 25.0e0
            c1min = log(h_charnock) / kappa
            do i = 1, maxiters
                f = c1_charnock - 2.0 * log(U10m)
                do j = 1, maxiters
                    c1 = (f + 2.0 * log(y1)) / kappa
                    ! looks like the check should use U10 instead of U
                    !   but note that a1 should be set = 0 in cycle beforehand
                    if(U <= 8.0e0) a1 = log(1.0 + c2_charnock * ((y1 / U10m)**3)) / kappa
                    c1 = max(c1 - a1, c1min)
                    y1 = c1
                end do
                z0_m = h_charnock * exp(-c1 * kappa)
                z0_m = max(z0_m,0.000015e0)
                U10m = U*alog(h_charnock/z0_m)/(alog(h/z0_m))
            end do

            ! --- define dynamic velocity in neutral conditions
            u_dyn0 = U10m / c1

        end subroutine
        ! --------------------------------------------------------------------------------


        SUBROUTINE  surf_flux(meteo, sfx, numerics, lu)
            type (sfxDataType), intent(out) :: sfx

            type (meteoDataType), intent(in) :: meteo
            type (numericsType), intent(in) :: numerics

            type (data_lutyp) lu

            real h              ! surface flux layer height [m]

            real z0_m, z0_t     ! aerodynamic & thermal roughness [m]
            real B              ! = ln(z0_m / z0_t) [n/d]
            real h0_m, h0_t     ! = h / z0_m, h / z0_h [n/d]

            real Re             ! roughness Reynolds number = u_dyn0 * z0 / nu [n/d]
            real u_dyn0         ! dynamic velocity in neutral conditions [m/s]

            real zeta           ! = z/L [n/d]
            real zeta_conv_lim  ! z/L critical value for matching free convection limit [n/d]

            real Rib            ! bulk Richardson number
            real Rib_conv_lim   ! Ri-bulk critical value for matching free convection limit [n/d]

            real psi_m, psi_h   ! universal functions [momentum] & [thermal]

            real U              ! wind velocity at h [m/s]
            real Tsemi          !


            integer lu_indx

            ! output ... !
            real an4, o, am
            real c4, c0, an0
            ! ...

            ! local unnamed !
            real x10, y10
            real al

            real q4, t4
            ! .....

            ! just local variables
            real f1, a1
            ! ....

            integer is_ocean
            integer i, j


            Tsemi = meteo%Tsemi

            U = meteo%U
            t4 = meteo%dT
            q4 = meteo%dQ
            h = meteo%h
            z0_m = meteo%z0_m

            if(z0_m < 0.0) then
                !> @brief Test - definition z0 of sea surface
                !> @details  if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
                is_ocean = 1

                call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)

                ! --- define relative height
                h0_m = h / z0_m
            else
                !     ......parameters from viscosity sublayer......
                is_ocean = 0

                ! --- define relative height
                h0_m = h / z0_m
                ! --- define dynamic velocity in neutral conditions
                u_dyn0 = U * kappa / log(h0_m)
            end if

            Re = u_dyn0 * z0_m / nu_air
            if(Re <= Re_rough_min) then
                B = B1_rough * alog(B3_rough * Re) + B2_rough
            else
                ! B4 takes into account Re value at z' ~ O(10) z0
                B = B4_rough * (Re**B2_rough)
            end if

            ! --- apply max restriction based on surface type
            if (is_ocean == 1) then
                B = min(B, B_max_ocean)
            else
                B = min(B, B_max_land)
            end if

            ! --- define roughness [thermal]
            z0_t = z0_m / exp(B)
            ! --- define relative height [thermal]
            h0_t = h / z0_t

            !     ......humidity stratification and ri-number......
            al = g / Tsemi
            Rib = al * h * (t4 + 0.61e0 * Tsemi * q4) / U**2

            ! --- define free convection transition zeta = z/L value
            call get_convection_lim(zeta_conv_lim, Rib_conv_lim, x10, y10, &
                h0_m, h0_t, B)

            if (Rib > 0.0e0) then
                !     ......stable stratification......
                !write (*,*) 'stable'

                Rib = min(Rib, Rib_max)
                call get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)

                f1 = beta_m * zeta
                am = 1.0 + f1
                o = 1.0/Pr_t_0_inv + f1

            else if (Rib < Rib_conv_lim) then
                !     ......strong instability.....
                !write (*,*) 'instability'

                call get_psi_convection(psi_m, psi_h, zeta, Rib, &
                        zeta_conv_lim, x10, y10, h0_m, h0_t, B, numerics%maxiters_convection)

                a1 = (zeta_conv_lim / zeta)**(1.0/3.0)
                am = a1 / y10
                o = a1 / (Pr_t_0_inv * x10)

            else if (Rib > -0.001e0) then

                !     ......nearly neutral......
                write (*,*) 'neutral'

                call get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)

                am = 1.0e0
                o = 1.0e0 / Pr_t_0_inv
            else
                !     ......week and semistrong instability......
                !write (*,*) 'semistrong'

                call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection)

                am = (1.0 - alpha_m * zeta)**(-0.25)
                o = 1.0/(Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta))
            end if

            !     ......computation of cu,co,k(h),alft
            c4=kappa/psi_m
            c0=kappa/psi_h
            an4=kappa*c4*U*h/am
            an0=am/o

            !     ......exit......
            sfx%zeta = zeta
            sfx%Rib = Rib
            sfx%Re = Re
            sfx%lnzuzt=B
            sfx%z0_m = z0_m
            sfx%z0_t = z0_t
            sfx%Rib_conv_lim = Rib_conv_lim
            sfx%cm=c4
            sfx%ch=c0
            sfx%ct=an4
            sfx%ckt=an0
            return

        END SUBROUTINE  surf_flux

    END MODULE drag3