    MODULE drag3
        USE param
        USE inputdata

        implicit none


        type, public:: data_in
	        real, public:: ws, dt, st, dq, cflh, z0in
        end type
	
        type, public:: data_outdef
	        real, public:: zl, ri, re, lnzuzt, zu, ztout, rith, cm, ch, ct, ckt
        end type
	
        type, public:: data_par
	        integer, public :: it=10
	    end type

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

    contains

        ! SUBROUTINE  surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, out, par1, lu1)
        SUBROUTINE surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, mas1_cflh, mas1_z0in, &
                masout_zl,  masout_ri, masout_re, masout_lnzuzt, masout_zu, masout_ztout,&
                masout_rith, masout_cm, masout_ch, masout_ct, masout_ckt,&
                par1, lu1,numst)

            real, dimension (numst) ::  mas1_w
            real, dimension (numst) ::  mas1_dt
            real, dimension (numst) ::  mas1_st
            real, dimension (numst) ::  mas1_dq
            real, dimension (numst) ::  mas1_cflh
            real, dimension (numst) ::  mas1_z0in

            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 (data_in)  in
            type (data_outdef) out
            type (data_par) par1
            type (data_lutyp) lu1
            integer i,numst

            do i=1,numst
                in%ws=mas1_w(i)
                in%dt=mas1_dt(i)
                in%st=mas1_st(i)
                in%dq=mas1_dq(i)
                in%cflh=mas1_cflh(i)
                in%z0in=mas1_z0in(i)

                !if ((in%ws == in%ws).and.(in%dt == in%dt).and.(in%st == in%st).and.(in%dq == in%dq)) then
                    CALL surf_flux(in, out, par1, lu1)
                !end if

                masout_zl(i)=out%zl
                masout_ri(i)=out%ri
                masout_re(i)=out%re
                masout_lnzuzt(i)=out%lnzuzt
                masout_zu(i)=out%zu
                masout_ztout(i)=out%ztout
                masout_rith(i)=out%rith
                masout_cm(i)=out%cm
                masout_ch(i)=out%ch
                masout_ct(i)=out%ct
                masout_ckt(i)=out%ckt
            enddo
        END SUBROUTINE  surf_fluxMAS

        SUBROUTINE  surf_flux(in, out, par, lu)
            type (data_in) , intent(in) :: in
            type (data_outdef) out
            type (data_par) par
            type (data_lutyp) lu

            real ws, dt, dq, cflh, z0in
	        integer it
            integer lu_indx
            real zl, ri, re, lnzuzt, zu, ztin, ri_th, cm, ch, ct, ckt
            real z0, d3, d0max, U10m, a1, y1, cimin, f, c1, u2, h0, u3, x7, d0, d00, zt, h00, ft0, an4, an5
            real al, Tsurf, Rib, q4, t4, u, r1, f0, f4, am, o, dd, x1, y0, x0, y10, a2ch, x10, p1, p0, h, d1, f1
            real d, c4, c1min, c0, c, b1, an0
            integer i, j

            ws=in%ws
            dt=in%dt
            Tsurf=in%st
            dq=in%dq
            cflh=in%cflh
            z0in=in%z0in
            it=par%it
	
            u=ws
            t4=dt
            q4=dq
            h=cflh
            z0=z0in

            AN5=(Pr_t_inf_inv / Pr_t_0_inv)**4
            D1=(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=(1.0E0-alpha_m*D1)**.25E0
            X10=(1.0E0-alpha_h*D1)**.5E0
            P1=2.0E0*ATAN(Y10)+ALOG((Y10-1.0E0)/(Y10+1.0E0))
            P0=ALOG((X10-1.0E0)/(X10+1.0E0))

            d3=0.0e0
            d0max=2.0e0

            if(z0 < 0.0e0) then
                !> @brief Test - definition z0 of sea surface
                !> @details  if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
                d0max = 8.0e0
                U10m=u
                a1=0.0e0
                y1=25.0e0
                c1min=alog(h1/1.0e0)/kappa
                do i=1,it
                    f=a2-2.0e0*alog(U10m)
                    do j=1,it
                        c1=(f+2.0e0*alog(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=alog(1.0e0+a3*((y1/U10m)**3))/kappa
                        c1=c1-a1
                        c1=max(c1,c1min)
                        y1=c1
                    end do
                    z0=h1*exp(-c1*kappa)
                    z0=max(z0,0.000015e0)
                    U10m=u*alog(h1/z0)/(alog(h/z0))
                end do
                h0=h/z0
                u3=U10m/c1
            else
                !     ......parameters from viscosity sublayer......
                h0=h/z0
                u3=u*kappa/alog(h0)
            end if

            x7=u3*z0/an
            if(x7 <= x8) then
                d0=an1*alog(al1*x7)+an2
            else
                d0=al2*(x7**0.45e0)
            end if
            !     ......humidity stratification and ri-number......
            al=g/Tsurf
            d0=min(d0,d0max)
            Rib=al*h*(t4+0.61e0*Tsurf*q4)/u**2
            d00=d0
            zt=z0/exp(d00)
            h00=h/zt
            ft0=alog(h00)
            !     ......definition of r-prim......
            an4=d1/h0
            an5=d1/h00
            if (abs(d0) < 1.0e-10) an5=an4
            an5=sqrt(1.0e0-g0*an5)
            an4=(1.0e0-alpha_m*an4)**0.25e0
            f0=alog((x10-1.0e0)*(an5+1.0e0)/((x10+1.0e0)*(an5-1.0e0)))/Pr_t_0_inv
            f4=2.0e0*(atan(y10)-atan(an4))+alog((y10-1.0e0)*(an4+1.0e0)/((y10+1.0e0)*(an4-1.0e0)))
            r1=d1*f0/(f4*f4)
            !     ......definition of dz,ta,fu,fo,fiu,fio......

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

                Rib=min(Rib,r0)
                f=alog(h0)
                f1=d0/f
                a1=b4*Rib
                a2ch=(f1+1.0e0)/Pr_t_0_inv-2.0e0*a1
                d3=f*(sqrt(a2ch**2+4.0e0*a1*(1.0e0-a1))-a2ch)/(2.0e0*b4*(1.0e0-a1))
                f1=b4*d3
                f4=f+f1
                f0=(f+d0)/Pr_t_0_inv+f1
                o=1.0e0/Pr_t_0_inv+f1
                am=1.0e0+f1
            else if (Rib < r1) then
                !     ......strong instability.....
                !write (*,*) 'instability'

                d3=d1
                do i=1,it+1
                    d=d3/h0
                    dd=d3/h00
                    if (abs(d0) < 1.0e-10) dd=d
                    a1=(d1/d3)**(1.0e0/3.0e0)
                    x0=sqrt(1.0e0-g0*dd)
                    y0=(1.0e0-alpha_m*d)**0.25e0
                    c=alog((x0+1.0e0)/(x0-1.0e0))
                    b1=-2.0e0*atan(y0)+alog((y0+1.0e0)/(y0-1.0e0))
                    f=3.0e0*(1.0e0-a1)
                    f4=f/y10+p1+b1
                    f0=(f/x10+p0+c)/Pr_t_0_inv
                    if (i == it+1) exit
                    d3=Rib*f4**2/f0
                end do
                am=a1/y10
                o=a1/(Pr_t_0_inv*x10)
            else if (Rib > -0.001e0) then

                !     ......nearly neutral......
                write (*,*) 'neutral'
                f4=alog(h0)
                f0=ft0/Pr_t_0_inv
                if (abs(d0) < 1.0e-10) f0=f4/Pr_t_0_inv
                am=1.0e0
                o=1.0e0/Pr_t_0_inv
            else
                !     ......week and semistrong instability......
                !write (*,*) 'semistrong'
                f1=alog(h0)
                if (abs(d0) < 1.0e-10) ft0=f1
                d3=Rib*Pr_t_0_inv*f1**2/ft0
                do i=1,it+1
                    d=d3/h0
                    dd=d3/h00
                    if (abs(d0) < 1.0e-10) dd=d
                    y1=(1.0e0-alpha_m*d3)**0.25e0
                    x1=sqrt(1.0e0-g0*d3)
                    y0=(1.0e0-alpha_m*d)**0.25e0
                    x0=sqrt(1.0e0-g0*dd)
                    y0=max(y0,1.000001e0)
                    x0=max(x0,1.000001e0)
                    f4=alog((y1-1.0e0)*(y0+1.0e0)/((y1+1.0e0)*(y0-1.0e0)))+2.0e0*(atan(y1)-atan(y0))
                    f0=alog((x1-1.0e0)*(x0+1.0e0)/((x1+1.0e0)*(x0-1.0e0)))/Pr_t_0_inv
                    if (i == it + 1) exit
                    d3=Rib*f4**2/f0
                end do
                am=(1.0e0-alpha_m*d3)**(-0.25e0)
                o=1.0e0/(Pr_t_0_inv*sqrt(1.0e0-g0*d3))
            end if

            !     ......computation of cu,co,k(h),alft
            c4=kappa/f4
            c0=kappa/f0
            an4=kappa*c4*u*h/am
            an0=am/o

            !     ......exit......
            out%zl=d3
            out%ri=Rib
            out%re=x7
            out%lnzuzt=d00
            out%zu=z0
            out%ztout=zt
            out%rith=r1
            out%cm=c4
            out%ch=c0
            out%ct=an4
            out%ckt=an0
            return

        END SUBROUTINE  surf_flux

    END MODULE drag3