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