Skip to content
Snippets Groups Projects
Commit a8269ccf authored by Evgeny Mortikov's avatar Evgeny Mortikov
Browse files

major code update to remove labels and goto statements

parent 7bd79379
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
compile2.sh
+ 4
− 4
View file @ a8269ccf
#!/bin/bash
rm drag_ddt.exe *.o
gfortran -c inputdata.f90
gfortran -c param.f90
gfortran -c drag3.f90
gfortran -c main_drag.f90
gfortran -c -Wuninitialized inputdata.f90
gfortran -c -Wuninitialized param.f90
gfortran -c -Wuninitialized drag3.f90
gfortran -c -Wuninitialized main_drag.f90
gfortran -o drag_ddt.exe main_drag.o drag3.o inputdata.o param.o
drag3.f90
+ 204
− 221
View file @ a8269ccf
MODULE drag3
USE param
USE inputdata
MODULE drag3
USE param
USE inputdata
implicit none
implicit none
type, public:: data_in
real, public:: ws, dt, st, dq, cflh, z0in
end type
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
end type
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
end type
type, public:: data_par
integer, public :: it=10
end type
type, public:: data_lutyp
integer, public :: lu_indx=1
end type
type, public:: data_lutyp
integer, public :: lu_indx=1
end type
contains
contains
! SUBROUTINE surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, out, par1, lu1)
SUBROUTINE surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, mas1_cflh, mas1_z0in, &
SUBROUTINE surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, mas1_cflh, mas1_z0in, &
masout_zl, masout_ri, masout_re, masout_lnzuzt, masout_zu, masout_ztout,&
masout_rith, masout_cm, masout_ch, masout_ct, masout_ckt,&
par1, lu1,numst)
......@@ -63,7 +62,9 @@
in%cflh=mas1_cflh(i)
in%z0in=mas1_z0in(i)
CALL surf_flux(in, out, par1, lu1)
!if ((in%ws == in%ws).and.(in%dt == in%dt).and.(in%st == in%st).and.(in%dq == in%dq)) then
CALL surf_flux(in, out, par1, lu1)
!end if
masout_zl(i)=out%zl
masout_ri(i)=out%ri
......@@ -80,209 +81,191 @@
END SUBROUTINE surf_fluxMAS
SUBROUTINE surf_flux(in, out, par, lu)
type (data_in) , intent(in) :: in
type (data_outdef) out
type (data_par) par
type (data_lutyp) lu
type (data_in) , intent(in) :: in
type (data_outdef) out
type (data_par) par
type (data_lutyp) lu
real ws, dt, dq, cflh, z0in
integer it
integer lu_indx
real zl, ri, re, lnzuzt, zu, ztin, ri_th, cm, ch, ct, ckt
real z0, d3, d0max, U10m, a1, y1, cimin, f, c1, u2, h0, u3, x7, d0, d00, zt, h00, ft0, an4, an5
real al, Tsurf, Rib, q4, t4, u, r1, f0, f4, am, o, dd, x1, y0, x0, y10, a2ch, x10, p1, p0, h, d1, f1
real d, c4, c1min, c0, c, b1, an0
integer i, j
real ws, dt, st, dq, cflh, z0in
integer it
integer lu_indx
real zl, ri, re, lnzuzt, zu, ztin, ri_th, cm, ch, ct, ckt
real z0, d3, d0max, u1, a1, y1, cimin, f, c1, u2, h0, u3, x7, d0, d00, zt, h00, ft0, an4, an5
real al, t1, r6, q4, t4, u, r1, f0, f4, 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
Tsurf=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
q4=dq
h=cflh
z0=z0in
AN5=(Pr_t_inf_inv / Pr_t_0_inv)**4
D1=(2.0E0*alpha_h-AN5*alpha_m-SQRT((AN5*alpha_m)**2+4.0E0*AN5*alpha_h*(alpha_h-alpha_m)))/(2.0E0*alpha_h**2)
Y10=(1.0E0-alpha_m*D1)**.25E0
X10=(1.0E0-alpha_h*D1)**.5E0
P1=2.0E0*ATAN(Y10)+ALOG((Y10-1.0E0)/(Y10+1.0E0))
P0=ALOG((X10-1.0E0)/(X10+1.0E0))
d3=0.0e0
d0max=2.0e0
if(z0 < 0.0e0) then
!> @brief Test - definition z0 of sea surface
!> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
d0max = 8.0e0
U10m=u
a1=0.0e0
y1=25.0e0
c1min=alog(h1/1.0e0)/kappa
do i=1,it
f=a2-2.0e0*alog(U10m)
do j=1,it
c1=(f+2.0e0*alog(y1))/kappa
! looks like the check should use U10 instead of U
! but note that a1 should be set = 0 in cycle beforehand
if(u <= 8.0e0) a1=alog(1.0e0+a3*((y1/U10m)**3))/kappa
c1=c1-a1
c1=max(c1,c1min)
y1=c1
end do
z0=h1*exp(-c1*kappa)
z0=max(z0,0.000015e0)
U10m=u*alog(h1/z0)/(alog(h/z0))
end do
h0=h/z0
u3=U10m/c1
else
! ......parameters from viscosity sublayer......
h0=h/z0
u3=u*kappa/alog(h0)
end if
x7=u3*z0/an
if(x7 <= x8) then
d0=an1*alog(al1*x7)+an2
else
d0=al2*(x7**0.45e0)
end if
! ......humidity stratification and ri-number......
al=g/Tsurf
d0=min(d0,d0max)
Rib=al*h*(t4+0.61e0*Tsurf*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
if (abs(d0) < 1.0e-10) an5=an4
an5=sqrt(1.0e0-g0*an5)
an4=(1.0e0-alpha_m*an4)**0.25e0
f0=alog((x10-1.0e0)*(an5+1.0e0)/((x10+1.0e0)*(an5-1.0e0)))/Pr_t_0_inv
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 (Rib > 0.0e0) then
! ......stable stratification......
!write (*,*) 'stable'
Rib=min(Rib,r0)
f=alog(h0)
f1=d0/f
a1=b4*Rib
a2ch=(f1+1.0e0)/Pr_t_0_inv-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)/Pr_t_0_inv+f1
o=1.0e0/Pr_t_0_inv+f1
am=1.0e0+f1
else if (Rib < r1) then
! ......strong instability.....
!write (*,*) 'instability'
d3=d1
do i=1,it+1
d=d3/h0
dd=d3/h00
if (abs(d0) < 1.0e-10) dd=d
a1=(d1/d3)**(1.0e0/3.0e0)
x0=sqrt(1.0e0-g0*dd)
y0=(1.0e0-alpha_m*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)/Pr_t_0_inv
if (i == it+1) exit
d3=Rib*f4**2/f0
end do
am=a1/y10
o=a1/(Pr_t_0_inv*x10)
else if (Rib > -0.001e0) then
! ......nearly neutral......
write (*,*) 'neutral'
f4=alog(h0)
f0=ft0/Pr_t_0_inv
if (abs(d0) < 1.0e-10) f0=f4/Pr_t_0_inv
am=1.0e0
o=1.0e0/Pr_t_0_inv
else
! ......week and semistrong instability......
!write (*,*) 'semistrong'
f1=alog(h0)
if (abs(d0) < 1.0e-10) ft0=f1
d3=Rib*Pr_t_0_inv*f1**2/ft0
do i=1,it+1
d=d3/h0
dd=d3/h00
if (abs(d0) < 1.0e-10) dd=d
y1=(1.0e0-alpha_m*d3)**0.25e0
x1=sqrt(1.0e0-g0*d3)
y0=(1.0e0-alpha_m*d)**0.25e0
x0=sqrt(1.0e0-g0*dd)
y0=max(y0,1.000001e0)
x0=max(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)))/Pr_t_0_inv
if (i == it + 1) exit
d3=Rib*f4**2/f0
end do
am=(1.0e0-alpha_m*d3)**(-0.25e0)
o=1.0e0/(Pr_t_0_inv*sqrt(1.0e0-g0*d3))
end if
! ......computation of cu,co,k(h),alft
c4=kappa/f4
c0=kappa/f0
an4=kappa*c4*u*h/am
an0=am/o
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))
d3=0.0e0
d0max=2.0e0
! ......exit......
out%zl=d3
out%ri=Rib
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
if(z0.lt.0.0e0) d0max=8.0e0
if(z0.lt.0.0e0) then
!> @brief Test - definition z0 of sea surface
!> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
END SUBROUTINE surf_flux
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
return
END SUBROUTINE surf_flux
END MODULE drag3
\ No newline at end of file
END MODULE drag3
\ No newline at end of file
main_drag.f90
+ 4
− 4
View file @ a8269ccf
......@@ -73,13 +73,13 @@
write(*,*) 'running code'
if (TEST==1) then
filename_in='MOSAiC.txt'
filename_in='data/MOSAiC.txt'
filename_out='out_MOSAiC.txt'
filename_in2='MOSAiC_zh.txt'
filename_in2='data/MOSAiC_zh.txt'
elseif (TEST==2) then
filename_in='Irgason1.txt'
filename_in='data/Irgason1.txt'
filename_out='out_IRGASON1.txt'
filename_in2='IRGASON_zh.txt'
filename_in2='data/IRGASON_zh.txt'
endif
open (1, file= filename_in, status ='old')
......
param.f90
+ 18
− 11
View file @ a8269ccf
module param
implicit real (a-h,o-z)
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
implicit none
! acceleration due to gravity [m/s^2]
real, parameter :: g = 9.81
! molecular Prandtl number (air)
real, parameter :: Pr_m = 0.71
! von Karman constant [n/d]
real, parameter :: kappa = 0.40
! inverse Prandtl (turbulent) number in neutral conditions
real, parameter :: Pr_t_0_inv = 1.15
! inverse Prandtl (turbulent) number in free convection
real, parameter :: Pr_t_inf_inv = 3.5
! stability function coeff. [= g4 & g10 in deprecated code]
real, parameter :: alpha_m = 16.0
real, parameter :: alpha_h = 16.0
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 :: al1=kappa*Pr_m
real, parameter :: g0=1.2
real, parameter :: al2=(.14e0*(30.0e0**an2))*(p4**.8e0)
real, parameter :: al2=(.14e0*(30.0e0**an2))*(Pr_m**.8e0)
real, parameter :: a2=alog(h1*(g/alfam))
real, parameter :: a3=betam*an*a2
real, parameter :: r0=.9e0/b4
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment