module DRAG3
   
      USE PRMT
      USE PARAM
      USE INPUTDATA
     
    !implicit real (A-H, O-Z)
    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
 

	contains
    
	subroutine  surf_flux(in, out, par)
    type (DATA_IN) , intent(in) :: in
    type (DATA_OUTDEF) out   
    type (DATA_PAR) par 

    real WS, DT, ST, DQ, CFLH, Z0IN
	integer IT
    real ZL, RI, RE, LNZUZT, ZU, ZTIN, RI_TH, CM, CH, CT, CKT
    real Z0, D3, D0MAX, U1, A1, Y1, CIMIN, H1, AP0, F, A2, C1, U2, H0, U3, X7, X8, AN1, AN2, D0, D00, ZT, H00, FT0, AN4, AN5
    real AL, AL2, AN, G, T1, R6, Q4, T4, U, G0, R1, F0, F4, A0, AM, O, DD, X1, Y0, X0, Z3, Y10, A2CH, X10, P1, P0, H, D1, F1
    real D, C4, C1MIN, C0, C, B1, AN0
    integer i, j, m
    WS=in%WS
    DT=in%DT
    ST=in%ST
    DQ=in%DQ
    CFLH=in%CFLH
    Z0IN=in%Z0IN	
    IT=par%IT
	
    U=WS
    T4=DT
    T1=ST
    Q4=DQ
    H=CFLH
    Z0=Z0IN
    
	  D3=0.0E0
      D0MAX=2.0E0
      !=DATA_IN%WS
      !4=DATA_IN%DT
      !4=DATA_IN%DQ
      !=DATA_IN%CFLH
      !0=DATA_IN%Z0
      IF(Z0.LT.0.0E0) D0MAX=8.0E0
      IF(Z0.LT.0.0E0) THEN
!     ......DEFINITION Z0 OF SEA SURFACE......
 !call Z0SEA module (RAMIL_DASHA)
                 !U1=U
      !A1=0.0E0
      !Y1=25.0E0
      !C1MIN=ALOG(H1/1.0E0)/AP0
      !DO 630 I=1,IT
      !F=A2-2.0E0*ALOG(U1)
      !DO 570 J=1,IT
      !C1=(F+2.0E0*ALOG(Y1))/AP0
      !IF(U.LE.8.0E0) A1=ALOG(1.0E0+A3*((Y1/U1)**3))/AP0
      !C1=C1-A1
      !C1=AMAX1(C1,C1MIN)
      !Y1=C1
  !570 CONTINUE
   !   Z0=H1*EXP(-C1*AP0)
   !   Z0=AMAX1(Z0,0.000015E0)
   !   U2=U*ALOG(H1/Z0)/(ALOG(H/Z0))
   !   U1=U2
  !630 CONTINUE
   !   J=1
   !   H0=H/Z0
   !   U3=U1/C1
    !  ELSE
!     ......PARAMETERS FROM VISCOSITY SUBLAYER......
    !  J=0
    !  H0=H/Z0
    !  U3=U*AP0/ALOG(H0)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      END IF
      X7=U3*Z0/AN
      IF(X7.LE.X8) THEN
      D0=AN1*ALOG(AL1*X7)+AN2
      ELSE
      D0=AL2*(X7**0.45E0)
      END IF
!     ......HUMIDITY STRATIFICATION AND RI-NUMBER......     
      ST=in%ST
      AL=G/T1
      D0=AMIN1(D0,D0MAX)
      R6=AL*H*(T4+0.61E0*T1*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
!C     IF(D0.EQ.0.0E0) AN5=AN4
      IF (ABS(D0).LT.1.0E-10) AN5=AN4
      AN5=SQRT(1.0E0-G0*AN5)
      AN4=(1.0E0-G4*AN4)**0.25E0
      F0=ALOG((X10-1.0E0)*(AN5+1.0E0)/((X10+1.0E0)*(AN5-1.0E0)))/A0
      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(R6.GT.0.0E0) GO TO 1460
      IF(R6.LT.R1) GO TO 1305
      IF(R6.GT.-0.001E0) THEN
!     ......NEARLY NEUTRAL......
      write (*,*) 'NEUTRAL'
      F4=ALOG(H0)
      F0=FT0/A0
!     IF(D0.EQ.0.0E0) F0=F4/A0
      IF (ABS(D0).LT.1.0E-10) F0=F4/A0
      AM=1.0E0
      O=1.0E0/A0
      GO TO 1570
      ELSE
!     ......WEEK AND SEMISTRONG INSTABILITY......
      write (*,*) 'SEMISTRONG'
      F1=ALOG(H0)
!     IF(D0.EQ.0.0E0)FT0=F1
      IF (ABS(D0).LT.1.0E-10) FT0=F1
      D3=R6*A0*F1**2/FT0
      M=1
 1245 DO 1300 I=1,IT
      D=D3/H0
      DD=D3/H00
!     IF(D0.EQ.0.0E0)DD=D
      IF (ABS(D0).LT.1.0E-10) DD=D
      Y1=(1.0E0-G4*D3)**0.25E0
      X1=SQRT(1.0E0-G0*D3)
      Y0=(1.0E0-G4*D)**0.25E0
      X0=SQRT(1.0E0-G0*DD)
      Y0=AMAX1(Y0,1.000001E0)
      X0=AMAX1(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)))/A0
      IF(M.NE.1) GO TO 1350
      Z3=R6*F4**2/F0
      D3=Z3
 1300 CONTINUE
      M=2
      GO TO 1245
 1350 AM=(1.0E0-G4*D3)**(-0.25E0)
      O=1.0E0/(A0*SQRT(1.0E0-G0*D3))
      GO TO 1570
      END IF
!     ......STRONG INSTABILITY.....
 1305 CONTINUE
      write (*,*) 'INSTABILITY'
      D3=D1
      M=1
 1355 DO 1410 I=1,IT
      D=D3/H0
      DD=D3/H00
!     IF(D0.EQ.0.0E0)DD=D
      IF (ABS(D0).LT.1.0E-10) DD=D
      A1=(D1/D3)**(1.0E0/3.0E0)
      X0=SQRT(1.0E0-G0*DD)
      Y0=(1.0E0-G4*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)/A0
      IF(M.NE.1) GO TO 1430
      Z3=R6*F4**2/F0
      D3=Z3
 1410 CONTINUE
      M=2
      GO TO 1355
 1430 AM=A1/Y10
      O=A1/(A0*X10)
      GO TO 1570
!     ......STABLE STRATIFICATION......
 1460 CONTINUE
      write (*,*) 'STABLE'
      R6=AMIN1(R6,R0)
      F=ALOG(H0)
      F1=D0/F
      A1=B4*R6
      A2CH=(F1+1.0E0)/A0-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)/A0+F1
      O=1.0E0/A0+F1
      AM=1.0E0+F1
 1570 CONTINUE
!     ......COMPUTATION OF CU,CO,K(H),ALFT
      C4=AP0/F4
      C0=AP0/F0
      AN4=AP0*C4*U*H/AM
      AN0=AM/O
!     ......EXIT......
  140 CONTINUE
      out%ZL=D3
      out%RI=R6
      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