MODULE drag3
USE param
USE inputdata
implicit none
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_par
integer, public :: it=10
end type
type, public:: data_lutyp
integer, public :: lu_indx=1
end type
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, &
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)
real, dimension (numst) :: mas1_w
real, dimension (numst) :: mas1_dt
real, dimension (numst) :: mas1_st
real, dimension (numst) :: mas1_dq
real, dimension (numst) :: mas1_cflh
real, dimension (numst) :: mas1_z0in
real, dimension (numst) :: masout_zl
real, dimension (numst) :: masout_ri
real, dimension (numst) :: masout_re
real, dimension (numst) :: masout_lnzuzt
real, dimension (numst) :: masout_zu
real, dimension (numst) :: masout_ztout
real, dimension (numst) :: masout_rith
real, dimension (numst) :: masout_cm
real, dimension (numst) :: masout_ch
real, dimension (numst) :: masout_ct
real, dimension (numst) :: masout_ckt
type (data_in) in
type (data_outdef) out
type (data_par) par1
type (data_lutyp) lu1
integer i,numst
do i=1,numst
in%ws=mas1_w(i)
in%dt=mas1_dt(i)
in%st=mas1_st(i)
in%dq=mas1_dq(i)
in%cflh=mas1_cflh(i)
in%z0in=mas1_z0in(i)
!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
masout_re(i)=out%re
masout_lnzuzt(i)=out%lnzuzt
masout_zu(i)=out%zu
masout_ztout(i)=out%ztout
masout_rith(i)=out%rith
masout_cm(i)=out%cm
masout_ch(i)=out%ch
masout_ct(i)=out%ct
masout_ckt(i)=out%ckt
enddo
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
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
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
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
! ......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
END SUBROUTINE surf_flux
END MODULE drag3