Skip to content
Snippets Groups Projects
Commit ba64fed7 authored by Victoria Suiazova's avatar Victoria Suiazova
Browse files

First commit

parent 8b8f21d1
Branches
Tags
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
COMP_DRAG.txt 0 → 100644
+ 7
− 0
View file @ ba64fed7
ifort -c inputdata.f90
ifort -c param.f90
ifort -c prmt.f90
ifort -c drag3.F
ifort -c DRAG.F
ifort -o drag.exe DRAG.o drag3.o inputdata.o param.o prmt.o
DRAG.F 0 → 100644
+ 33
− 0
View file @ ba64fed7
PROGRAM DRAG
USE INPUTDATA
USE PARAM
open (1, file= '2016_inp.txt', status ='old')
open (2, file='2016_outDRAG.txt', status='new')
do i=1,n
read (1,*) ws, deltaT, semisumT
AR1(1) = ws
AR1(2) = deltaT
AR1(3) = semisumT
AR1(4) = 0
AR1(5) = 9
AR1(6) = 0.02
AN5=(A6/A0)**4
D1=(2.0E0*G10-AN5*G4-SQRT((AN5*G4)**2+4.0E0*AN5*G10*(G10-G4)))
*/(2.0E0*G10**2)
Y10=(1.0E0-G4*D1)**.25E0
X10=(1.0E0-G10*D1)**.5E0
P1=2.0E0*ATAN(Y10)+ALOG((Y10-1.0E0)/(Y10+1.0E0))
P0=ALOG((X10-1.0E0)/(X10+1.0E0))
CALL DRAG3
HFX=-RO*CP*U*T4*C0
MFX=-RO*C4*U*U
betta=G/T1
zL=(betta*AKA*H*C0*T4)/(C4*C4*U*U)
write (2,*) i, HFX, MFX, ZL
write (*,*) 'happy end'
enddo
END PROGRAM
data/2016_inp.txt 0 → 100644
+ 7094
− 0
View file @ ba64fed7
This diff is collapsed.
drag3.F 0 → 100644
+ 204
− 0
View file @ ba64fed7
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
inputdata.f90 0 → 100644
+ 27
− 0
View file @ ba64fed7
module INPUTDATA
REAl, DIMENSION (6) :: AR1
REAl, DIMENSION (11) :: AR2
INTEGER, PARAMETER :: IT =1
INTEGER, PARAMETER :: N=7094
REAL HFX, MFX, zL, betta
REAL U, T4,C0,C4, T1,H
REAL ws, deltaT, semisumT
!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*====================================================================
end module INPUTDATA
\ No newline at end of file
param.f90 0 → 100644
+ 50
− 0
View file @ ba64fed7
module PARAM
IMPLICIT REAL (A-H,O-Z)
real, parameter :: RO = 1.2
real, parameter :: CP =8.4
real, parameter :: AKA=.40E0
real, parameter :: AP0=.40E0
real, parameter :: G=9.81E0
real, parameter :: A0=1.15E0
real, parameter :: A6=3.5E0
real, parameter :: G4=16.0E0
real, parameter :: G10=16.0E0
real, parameter :: B4=4.7E0
real, parameter :: ALFAM=.0144E0
real, parameter :: BETAM=.111E0
real, parameter :: AN=.000015E0
real, parameter :: P4=.71E0
real, parameter :: H1=10.0E0
real, parameter :: X8=16.3E0
real, parameter :: AN1=5.0E0/6.0E0
real, parameter :: AN2=.45E0
real, parameter :: AL1=AKA*P4
real, parameter :: G0=1.2
real, parameter :: AL2=(.14E0*(30.0E0**AN2))*(P4**.8E0)
real, parameter :: A2=ALOG(H1*(G/ALFAM))
real, parameter :: A3=BETAM*AN*A2
real, parameter :: R0=.9E0/B4
!C* real, parameter :: AN5=(A6/A0)**4
!C* real, parameter :: D1=(2.0E0*G10-AN5*G4-SQRT((AN5*G4)**2+4.0E0*AN5*G10*(G10-G4)))/(2.0E0*G10**2)
!C* real, parameter :: Y10=(1.0E0-G4*D1)**.25E0
!C* real, parameter :: X10=(1.0E0-G10*D1)**.5E0
!C* real, parameter :: P1=2.0E0*ATAN(Y10)+ALOG((Y10-1.0E0)/(Y10+1.0E0))
!C* real, parameter :: P0=ALOG((X10-1.0E0)/(X10+1.0E0))
! AN5=(A6/A0)**4
! D1=(2.0E0*G10-AN5*G4-SQRT((AN5*G4)**2+4.0E0*AN5*G10*(G10-G4)))/(2.0E0*G10**2)
! Y10=(1.0E0-G4*D1)**.25E0
! X10=(1.0E0-G10*D1)**.5E0
! P1=2.0E0*ATAN(Y10)+ALOG((Y10-1.0E0)/(Y10+1.0E0))
! P0=ALOG((X10-1.0E0)/(X10+1.0E0))
end module PARAM
prmt.f90 0 → 100644
+ 39
− 0
View file @ ba64fed7
module PRMT
implicit none
integer, parameter :: NLON = 180
integer, parameter :: NLAT = 121
integer, parameter :: NLEV = 21
integer, parameter :: NLATM1 = NLAT - 1
integer, parameter :: NLONP1 = NLON + 1
integer, parameter :: NLONP2 = NLON + 2
integer, parameter :: NLEVM1 = NLEV - 1
integer, parameter :: LBUF = NLONP2 * NLEV * 5
integer, parameter :: LENTX = NLONP2 * NLAT * NLEV * 10 + LBUF * 2
integer, parameter :: NR21 = NLEV + 1
integer, parameter :: KL = NLEV
integer, parameter :: KST = 17
integer, parameter :: KLM = 21
integer, parameter :: NTR = 10
integer, parameter :: NTRY = 3
integer, parameter :: NFFT = NLON+MOD(NLON,32)+32
integer, parameter :: LENTXX = NLONP2 * NLAT * NLEV * NTR * 2 + NLONP2 * NLEV * NTR * 2
integer, parameter :: IA = 1
integer, parameter :: LakeModel = 0 ! 0 - off the Lake model, old lakes parametrization with 13 surfaces tupes from VEG
! 1 - on the Lake model, online coupling the Lake model with 14 surfaces types from VEGNEW and
! lakes depths from ADLake
! 2 - on the Lake model, old lakes parametrization and offline coupling the Lake model
! with 14 surfaces types from VEGNEW and a depths from ADLake
end module PRMT
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment