Skip to content
Snippets Groups Projects
Commit 5215915f authored by Бычкова Виктория's avatar Бычкова Виктория
Browse files

low regidtr

parent 84976472
No related branches found
No related tags found
No related merge requests found
...@@ -6,5 +6,5 @@ gfortran -c param.f90 ...@@ -6,5 +6,5 @@ gfortran -c param.f90
gfortran -c prmt.f90 gfortran -c prmt.f90
gfortran -c drag3.f90 gfortran -c drag3.f90
gfortran -c main_drag.f90 gfortran -c main_drag.f90
gfortran -o drag_ddt.exe main_drag.o drag3.o inputdata.o param.o prmt.o gfortran -o drag_ddt.exe main_drag.o drag3.o inputdata.o param.o
module DRAG3 MODULE drag3
USE PRMT USE param
USE PARAM USE inputdata
USE INPUTDATA
!implicit real (A-H, O-Z) !implicit real (a-h, o-z)
implicit none implicit none
type, public:: DATA_IN type, public:: data_in
real, public:: WS, DT, ST, DQ, CFLH, Z0IN real, public:: ws, dt, st, dq, cflh, z0in
end type end type
type, public:: DATA_OUTDEF type, public:: data_outdef
real, public:: ZL, RI, RE, LNZUZT, ZU, ZTOUT, RITH, CM, CH, CT, CKT real, public:: zl, ri, re, lnzuzt, zu, ztout, rith, cm, ch, ct, ckt
end type end type
type, public:: DATA_PAR type, public:: data_par
integer, public :: IT=10 integer, public :: it=10
end type end type
contains contains
subroutine surf_flux(in, out, par) SUBROUTINE surf_flux(in, out, par)
type (DATA_IN) , intent(in) :: in type (data_in) , intent(in) :: in
type (DATA_OUTDEF) out type (data_outdef) out
type (DATA_PAR) par type (data_par) par
real WS, DT, ST, DQ, CFLH, Z0IN real ws, dt, st, dq, cflh, z0in
integer IT integer it
real ZL, RI, RE, LNZUZT, ZU, ZTIN, RI_TH, CM, CH, CT, CKT 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 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 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 real d, c4, c1min, c0, c, b1, an0
integer i, j, m integer i, j, m
WS=in%WS ws=in%ws
DT=in%DT dt=in%dt
ST=in%ST st=in%st
DQ=in%DQ dq=in%dq
CFLH=in%CFLH cflh=in%cflh
Z0IN=in%Z0IN z0in=in%z0in
IT=par%IT it=par%it
U=WS u=ws
T4=DT t4=dt
T1=ST t1=st
Q4=DQ q4=dq
H=CFLH h=cflh
Z0=Z0IN z0=z0in
D3=0.0E0 d3=0.0e0
D0MAX=2.0E0 d0max=2.0e0
!=DATA_IN%WS !=data_in%ws
!4=DATA_IN%DT !4=data_in%dt
!4=DATA_IN%DQ !4=data_in%dq
!=DATA_IN%CFLH !=data_in%cflh
!0=DATA_IN%Z0 !0=data_in%z0
IF(Z0.LT.0.0E0) D0MAX=8.0E0 if(z0.lt.0.0e0) d0max=8.0e0
IF(Z0.LT.0.0E0) THEN if(z0.lt.0.0e0) then
! ......DEFINITION Z0 OF SEA SURFACE...... ! ......definition z0 of sea surface......
!call Z0SEA module (RAMIL_DASHA) !call z0sea module (ramil_dasha)
!U1=U !u1=u
!A1=0.0E0 !a1=0.0e0
!Y1=25.0E0 !y1=25.0e0
!C1MIN=ALOG(H1/1.0E0)/AP0 !c1min=alog(h1/1.0e0)/ap0
!DO 630 I=1,IT !do 630 i=1,it
!F=A2-2.0E0*ALOG(U1) !f=a2-2.0e0*alog(u1)
!DO 570 J=1,IT !do 570 j=1,it
!C1=(F+2.0E0*ALOG(Y1))/AP0 !c1=(f+2.0e0*alog(y1))/ap0
!IF(U.LE.8.0E0) A1=ALOG(1.0E0+A3*((Y1/U1)**3))/AP0 !if(u.le.8.0e0) a1=alog(1.0e0+a3*((y1/u1)**3))/ap0
!C1=C1-A1 !c1=c1-a1
!C1=AMAX1(C1,C1MIN) !c1=amax1(c1,c1min)
!Y1=C1 !y1=c1
!570 CONTINUE !570 continue
! Z0=H1*EXP(-C1*AP0) ! z0=h1*exp(-c1*ap0)
! Z0=AMAX1(Z0,0.000015E0) ! z0=amax1(z0,0.000015e0)
! U2=U*ALOG(H1/Z0)/(ALOG(H/Z0)) ! u2=u*alog(h1/z0)/(alog(h/z0))
! U1=U2 ! u1=u2
!630 CONTINUE !630 continue
! J=1 ! j=1
! H0=H/Z0 ! h0=h/z0
! U3=U1/C1 ! u3=u1/c1
! ELSE ! else
! ......PARAMETERS FROM VISCOSITY SUBLAYER...... ! ......parameters from viscosity sublayer......
! J=0 ! j=0
! H0=H/Z0 ! h0=h/z0
! U3=U*AP0/ALOG(H0) ! u3=u*ap0/alog(h0)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END IF end if
X7=U3*Z0/AN x7=u3*z0/an
IF(X7.LE.X8) THEN if(x7.le.x8) then
D0=AN1*ALOG(AL1*X7)+AN2 d0=an1*alog(al1*x7)+an2
ELSE else
D0=AL2*(X7**0.45E0) d0=al2*(x7**0.45e0)
END IF end if
! ......HUMIDITY STRATIFICATION AND RI-NUMBER...... ! ......humidity stratification and ri-number......
ST=in%ST st=in%st
AL=G/T1 al=g/t1
D0=AMIN1(D0,D0MAX) d0=amin1(d0,d0max)
R6=AL*H*(T4+0.61E0*T1*Q4)/U**2 r6=al*h*(t4+0.61e0*t1*q4)/u**2
D00=D0 d00=d0
ZT=Z0/EXP(D00) zt=z0/exp(d00)
H00=H/ZT h00=h/zt
FT0=ALOG(H00) ft0=alog(h00)
! ......DEFINITION OF R-PRIM...... ! ......definition of r-prim......
AN4=D1/H0 an4=d1/h0
AN5=D1/H00 an5=d1/h00
!C IF(D0.EQ.0.0E0) AN5=AN4 !c if(d0.eq.0.0e0) an5=an4
IF (ABS(D0).LT.1.0E-10) AN5=AN4 if (abs(d0).lt.1.0e-10) an5=an4
AN5=SQRT(1.0E0-G0*AN5) an5=sqrt(1.0e0-g0*an5)
AN4=(1.0E0-G4*AN4)**0.25E0 an4=(1.0e0-g4*an4)**0.25e0
F0=ALOG((X10-1.0E0)*(AN5+1.0E0)/((X10+1.0E0)*(AN5-1.0E0)))/A0 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))) 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) r1=d1*f0/(f4*f4)
! ......DEFINITION OF DZ,TA,FU,FO,FIU,FIO...... ! ......definition of dz,ta,fu,fo,fiu,fio......
IF(R6.GT.0.0E0) GO TO 1460 if(r6.gt.0.0e0) go to 1460
IF(R6.LT.R1) GO TO 1305 if(r6.lt.r1) go to 1305
IF(R6.GT.-0.001E0) THEN if(r6.gt.-0.001e0) then
! ......NEARLY NEUTRAL...... ! ......nearly neutral......
write (*,*) 'NEUTRAL' write (*,*) 'neutral'
F4=ALOG(H0) f4=alog(h0)
F0=FT0/A0 f0=ft0/a0
! IF(D0.EQ.0.0E0) F0=F4/A0 ! if(d0.eq.0.0e0) f0=f4/a0
IF (ABS(D0).LT.1.0E-10) F0=F4/A0 if (abs(d0).lt.1.0e-10) f0=f4/a0
AM=1.0E0 am=1.0e0
O=1.0E0/A0 o=1.0e0/a0
GO TO 1570 go to 1570
ELSE else
! ......WEEK AND SEMISTRONG INSTABILITY...... ! ......week and semistrong instability......
write (*,*) 'SEMISTRONG' write (*,*) 'semistrong'
F1=ALOG(H0) f1=alog(h0)
! IF(D0.EQ.0.0E0)FT0=F1 ! if(d0.eq.0.0e0)ft0=f1
IF (ABS(D0).LT.1.0E-10) FT0=F1 if (abs(d0).lt.1.0e-10) ft0=f1
D3=R6*A0*F1**2/FT0 d3=r6*a0*f1**2/ft0
M=1 m=1
1245 DO 1300 I=1,IT 1245 do 1300 i=1,it
D=D3/H0 d=d3/h0
DD=D3/H00 dd=d3/h00
! IF(D0.EQ.0.0E0)DD=D ! if(d0.eq.0.0e0)dd=d
IF (ABS(D0).LT.1.0E-10) DD=D if (abs(d0).lt.1.0e-10) dd=d
Y1=(1.0E0-G4*D3)**0.25E0 y1=(1.0e0-g4*d3)**0.25e0
X1=SQRT(1.0E0-G0*D3) x1=sqrt(1.0e0-g0*d3)
Y0=(1.0E0-G4*D)**0.25E0 y0=(1.0e0-g4*d)**0.25e0
X0=SQRT(1.0E0-G0*DD) x0=sqrt(1.0e0-g0*dd)
Y0=AMAX1(Y0,1.000001E0) y0=amax1(y0,1.000001e0)
X0=AMAX1(X0,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)) 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 f0=alog((x1-1.0e0)*(x0+1.0e0)/((x1+1.0e0)*(x0-1.0e0)))/a0
IF(M.NE.1) GO TO 1350 if(m.ne.1) go to 1350
Z3=R6*F4**2/F0 z3=r6*f4**2/f0
D3=Z3 d3=z3
1300 CONTINUE 1300 continue
M=2 m=2
GO TO 1245 go to 1245
1350 AM=(1.0E0-G4*D3)**(-0.25E0) 1350 am=(1.0e0-g4*d3)**(-0.25e0)
O=1.0E0/(A0*SQRT(1.0E0-G0*D3)) o=1.0e0/(a0*sqrt(1.0e0-g0*d3))
GO TO 1570 go to 1570
END IF end if
! ......STRONG INSTABILITY..... ! ......strong instability.....
1305 CONTINUE 1305 continue
write (*,*) 'INSTABILITY' write (*,*) 'instability'
D3=D1 d3=d1
M=1 m=1
1355 DO 1410 I=1,IT 1355 do 1410 i=1,it
D=D3/H0 d=d3/h0
DD=D3/H00 dd=d3/h00
! IF(D0.EQ.0.0E0)DD=D ! if(d0.eq.0.0e0)dd=d
IF (ABS(D0).LT.1.0E-10) DD=D if (abs(d0).lt.1.0e-10) dd=d
A1=(D1/D3)**(1.0E0/3.0E0) a1=(d1/d3)**(1.0e0/3.0e0)
X0=SQRT(1.0E0-G0*DD) x0=sqrt(1.0e0-g0*dd)
Y0=(1.0E0-G4*D)**0.25E0 y0=(1.0e0-g4*d)**0.25e0
C=ALOG((X0+1.0E0)/(X0-1.0E0)) c=alog((x0+1.0e0)/(x0-1.0e0))
B1=-2.0E0*ATAN(Y0)+ALOG((Y0+1.0E0)/(Y0-1.0E0)) b1=-2.0e0*atan(y0)+alog((y0+1.0e0)/(y0-1.0e0))
F=3.0E0*(1.0E0-A1) f=3.0e0*(1.0e0-a1)
F4=F/Y10+P1+B1 f4=f/y10+p1+b1
F0=(F/X10+P0+C)/A0 f0=(f/x10+p0+c)/a0
IF(M.NE.1) GO TO 1430 if(m.ne.1) go to 1430
Z3=R6*F4**2/F0 z3=r6*f4**2/f0
D3=Z3 d3=z3
1410 CONTINUE 1410 continue
M=2 m=2
GO TO 1355 go to 1355
1430 AM=A1/Y10 1430 am=a1/y10
O=A1/(A0*X10) o=a1/(a0*x10)
GO TO 1570 go to 1570
! ......STABLE STRATIFICATION...... ! ......stable stratification......
1460 CONTINUE 1460 continue
write (*,*) 'STABLE' write (*,*) 'stable'
R6=AMIN1(R6,R0) r6=amin1(r6,r0)
F=ALOG(H0) f=alog(h0)
F1=D0/F f1=d0/f
A1=B4*R6 a1=b4*r6
A2CH=(F1+1.0E0)/A0-2.0E0*A1 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)) d3=f*(sqrt(a2ch**2+4.0e0*a1*(1.0e0-a1))-a2ch)/(2.0e0*b4*(1.0e0-a1))
F1=B4*D3 f1=b4*d3
F4=F+F1 f4=f+f1
F0=(F+D0)/A0+F1 f0=(f+d0)/a0+f1
O=1.0E0/A0+F1 o=1.0e0/a0+f1
AM=1.0E0+F1 am=1.0e0+f1
1570 CONTINUE 1570 continue
! ......COMPUTATION OF CU,CO,K(H),ALFT ! ......computation of cu,co,k(h),alft
C4=AP0/F4 c4=ap0/f4
C0=AP0/F0 c0=ap0/f0
AN4=AP0*C4*U*H/AM an4=ap0*c4*u*h/am
AN0=AM/O an0=am/o
! ......EXIT...... ! ......exit......
140 CONTINUE 140 continue
out%ZL=D3 out%zl=d3
out%RI=R6 out%ri=r6
out%RE=X7 out%re=x7
out%LNZUZT=D00 out%lnzuzt=d00
out%ZU=Z0 out%zu=z0
out%ZTOUT=ZT out%ztout=zt
out%RITH=R1 out%rith=r1
out%CM=C4 out%cm=c4
out%CH=C0 out%ch=c0
out%CT=AN4 out%ct=an4
out%CKT=AN0 out%ckt=an0
return return
end subroutine surf_flux END SUBROUTINE surf_flux
end module DRAG3 END MODULE drag3
\ No newline at end of file \ No newline at end of file
program main_ddt PROGRAM main_ddt
USE PRMT
USE PARAM
USE INPUTDATA
use DRAG3
type (DATA_IN):: DATA_IN1 USE param
USE inputdata
USE drag3
type (DATA_OUTDEF) :: DATA_OUTDEF1 type (data_in):: data_in1
type (data_outdef) :: data_outdef1
type (DATA_PAR) :: DATA_PAR1
type (data_par) :: data_par1
OPEN (1,FILE='4_ddt.txt')
DO I = 1,1000000 open (1,file='4_ddt.txt')
READ (1,*,END=100) DATA_IN1%WS, DATA_IN1%DT
do i = 1,1000000
read (1,*,end=100) data_in1%ws, data_in1%dt
call surf_flux(DATA_IN1, DATA_OUTDEF1, DATA_PAR1)
CALL surf_flux(data_in1, data_outdef1, data_par1)
ENDDO
100 CONTINUE enddo
10 FORMAT (4I4,5F7.1,F7.4,F7.1) 100 continue
20 FORMAT (4I4,5F7.1,F7.4,F7.1)
10 format (4i4,5f7.1,f7.4,f7.1)
20 format (4i4,5f7.1,f7.4,f7.1)
stop stop
end program END PROGRAM
\ No newline at end of file
prmt.mod 0 → 100644
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment