      SUBROUTINE DRAG3
      USE PRMT
      USE PARAM
      USE INPUTDATA
      IMPLICIT REAL (A-H,O-Z)
C
C*====================================================================
C*     .....DEFENITION OF DRAG AND HEAT EXCHANGE COEFFICIENTS......  =
C*     DETAILS OF ALGORITM ARE GIVEN IN:                             =
C*     A.L.KAZAKOV,V.N.LYKOSSOV,"TRUDY ZAP.SIB.NII",1982,N.55,3-20   =
C*                    INPUT DATA:
C*     AR1(1) - ABS(WIND VELOCITY) AT CONSTANT FLUX LAYER            =
C*                                 (CFL) HIGHT (M/S)                 =
C*     AR1(2) - DIFFERENCE BETWEEN POTENTIAL TEMPERATURE AT CFL HIGHT=
C*                             AND AT SURFACE  ( DEG. K)             =
C*     AR1(3) - SEMI-SUM OF POTENTIAL TEMPERATURE AT CFL HIGHT AND   =
C*                             AND AT SURFACE  ( DEG. K)             =
C*     AR1(4) - DIFFERENCE BETWEEN HUMIDITY AT CFL HIGHT             =
C*                             AND A SURFACE   ( GR/GR )             =
C*     AR1(5) - CFL HIGHT ( M )                                      =
C*     AR1(6) - ROUGHNESS OF SURFACE ( M ); FOR SEA SURFACE PUT -1   =
C*      IT    - NUMBER OF ITERATIONS                                 =
C*                     OUTPUT DATA:
C*     AR2(1) - NON-DIMENSIONAL CFL HIGHT                            =
C*     AR2(2) - RICHARDSON NUMBER                                    =
C*     AR2(3) - REYNODS NUMBER                                       =
C*     AR2(4) - LN(ZU/ZT)                                            =
C*     AR2(5) - DYNAMICAL ROUGHNESS ZU (M)                           =
C*     AR2(6) - THERMAL   ROUGHNESS ZT (M)                           =
C*     AR2(7) - CRITICAL RICHARDSON NUMBER                           =
C*     AR2(8) - TRANSFER COEFFICIENT FOR MOMENTUM                    =
C*     AR2(9) - TRANSFER COEFFICIENT FR HEAT                         =
C*     AR2(10)- COEFFICIENT OF TURBULENCE (KM) AT CFL HIGHT (M**2/S) =
C*     AR2(11)- ALFT=KT/KM ( KT-COEFFICIENT OF TURBULENCE FOR HEAT)  =
C*     COMMENT: DRAG COEFFICIENT =          AR2(8)*AR2(8)            =
C*              HEAT EXCHANGE COEFFICIENT = AR2(8)*AR2(9)            =
C*====================================================================
      D3=0.0E0
      D0MAX=2.0E0
      U=AR1(1)
      T4=AR1(2)
      Q4=AR1(4)
      H=AR1(5)
      Z0=AR1(6)
      IF(Z0.LT.0.0E0) D0MAX=8.0E0
      IF(Z0.LT.0.0E0) THEN
C*     ......DEFINITION Z0 OF SEA SURFACE......
      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
C*     ......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
C*     ......HUMIDITY STRATIFICATION AND RI-NUMBER......
      T1=AR1(3)
      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)
C*     ......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)
C*     ......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
C*     ......NEARLY NEUTRAL......
      F4=ALOG(H0)
      F0=FT0/A0
C     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
C*     ......WEEK AND SEMISTRONG INSTABILITY......
      F1=ALOG(H0)
C     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
C     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
C*     ......STRONG INSTABILITY.....
 1305 CONTINUE
      D3=D1
      M=1
 1355 DO 1410 I=1,IT
      D=D3/H0
      DD=D3/H00
C     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
C*     ......STABLE STRATIFICATION......
 1460 CONTINUE
      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
C*     ......COMPUTATION OF CU,CO,K(H),ALFT
      C4=AP0/F4
      C0=AP0/F0
      AN4=AP0*C4*U*H/AM
      AN0=AM/O
C*     ......EXIT......
  140 CONTINUE
      AR2(1)=D3
      AR2(2)=R6
      AR2(3)=X7
      AR2(4)=D00
      AR2(5)=Z0
      AR2(6)=ZT
      AR2(7)=R1
      AR2(8)=C4
      AR2(9)=C0
      AR2(10)=AN4
      AR2(11)=AN0
      RETURN
      END SUBROUTINE
