diff --git a/compile2.sh b/compile2.sh
index b6966b86a84787045b697ffc2d8621c5353e8a8b..e1433f88e7ddfcdc6627204ad266a985b7e85a89 100644
--- a/compile2.sh
+++ b/compile2.sh
@@ -6,5 +6,5 @@ gfortran -c param.f90
gfortran -c prmt.f90
gfortran -c drag3.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
diff --git a/drag3.f90 b/drag3.f90
index 860abae30aa7368ccfe81c9fab1175f0c27fa156..039a999c735642db86d2f47cbcd2b1711d7b5715 100644
--- a/drag3.f90
+++ b/drag3.f90
@@ -1,224 +1,223 @@
- module DRAG3
+ MODULE drag3
- USE PRMT
- USE PARAM
- USE INPUTDATA
+ USE param
+ USE inputdata
- !implicit real (A-H, O-Z)
+ !implicit real (a-h, o-z)
implicit none
- type, public:: DATA_IN
- real, public:: WS, DT, ST, DQ, CFLH, Z0IN
+ 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
+ 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
+ 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
+ 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
+ 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
+ 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
+ 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)
+ 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
+ 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
\ No newline at end of file
+ END SUBROUTINE surf_flux
+ END MODULE drag3
\ No newline at end of file
diff --git a/main_drag.f90 b/main_drag.f90
index 4e50962853889feedb19161265d4b49f9526bebb..6c809aff54932813bc9aa783c8e3c3ac0d0358bf 100644
--- a/main_drag.f90
+++ b/main_drag.f90
@@ -1,33 +1,33 @@
-program main_ddt
- USE PRMT
- USE PARAM
- USE INPUTDATA
- use DRAG3
+ PROGRAM main_ddt
+
+ USE param
+ USE inputdata
+ USE drag3
- type (DATA_IN):: DATA_IN1
+ type (data_in):: data_in1
- type (DATA_OUTDEF) :: DATA_OUTDEF1
+ type (data_outdef) :: data_outdef1
- type (DATA_PAR) :: DATA_PAR1
+ type (data_par) :: data_par1
- OPEN (1,FILE='4_ddt.txt')
+ open (1,file='4_ddt.txt')
- DO I = 1,1000000
- 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
+ enddo
-100 CONTINUE
+100 continue
-10 FORMAT (4I4,5F7.1,F7.4,F7.1)
-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
-end program
+END PROGRAM
\ No newline at end of file
diff --git a/prmt.mod b/prmt.mod
new file mode 100644
index 0000000000000000000000000000000000000000..fbc5c106bbef1e8d13147a0421ec45fb0d9c8dd0
Binary files /dev/null and b/prmt.mod differ