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

major code update

parent a8269ccf
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
drag3.f90
+ 429
− 195
View file @ c1e705f6
MODULE drag3 MODULE drag3
USE param USE param
USE inputdata
implicit none implicit none
type, public :: meteoDataType
real :: h ! constant flux layer height [m]
real :: U ! abs(wind speed) at 'h' [m/s]
real :: dT ! difference between potential temperature at 'h' and at surface [K]
real :: Tsemi ! semi-sum of potential temperature at 'h' and at surface [K]
real :: dQ ! difference between humidity at 'h' and at surface [g/g]
real :: z0_m ! surface aerodynamic roughness (should be < 0 for water bodies surface)
end type
type, public :: meteoDataVecType
real, allocatable :: h(:) ! constant flux layer height [m]
real, allocatable :: U(:) ! abs(wind speed) at 'h' [m/s]
real, allocatable :: dT(:) ! difference between potential temperature at 'h' and at surface [K]
real, allocatable :: Tsemi(:) ! semi-sum of potential temperature at 'h' and at surface [K]
real, allocatable :: dQ(:) ! difference between humidity at 'h' and at surface [g/g]
real, allocatable :: z0_m(:) ! surface aerodynamic roughness (should be < 0 for water bodies surface)
end type
type, public:: data_in type, public :: sfxDataType
real, public:: ws, dt, st, dq, cflh, z0in real :: zeta ! = z/L [n/d]
real :: Rib ! bulk Richardson number [n/d]
real :: Re ! Reynolds number [n/d]
real :: lnzuzt
real :: z0_m ! aerodynamic roughness [m]
real :: z0_t ! thermal roughness [m]
real :: Rib_conv_lim ! Ri-bulk convection critical value [n/d]
real :: cm
real :: ch
real :: ct
real :: ckt
end type end type
type, public:: data_outdef type, public :: numericsType
real, public:: zl, ri, re, lnzuzt, zu, ztout, rith, cm, ch, ct, ckt integer :: maxiters_convection = 10 ! maximum (actual) number of iterations in convection
integer :: maxiters_charnock = 10 ! maximum (actual) number of iterations in charnock roughness
end type end type
type, public:: data_par
integer, public :: it=10
end type
type, public:: data_lutyp type, public:: data_lutyp
integer, public :: lu_indx=1 integer, public :: lu_indx=1
...@@ -23,18 +46,13 @@ ...@@ -23,18 +46,13 @@
contains contains
! SUBROUTINE surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, out, par1, lu1) subroutine surf_flux_vec(meteo, &
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_zl, masout_ri, masout_re, masout_lnzuzt, masout_zu, masout_ztout,&
masout_rith, masout_cm, masout_ch, masout_ct, masout_ckt,& masout_rith, masout_cm, masout_ch, masout_ct, masout_ckt,&
par1, lu1,numst) numerics, lu1,numst)
real, dimension (numst) :: mas1_w type (meteoDataVecType), intent(in) :: meteo
real, dimension (numst) :: mas1_dt type (numericsType), intent(in) :: numerics
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_zl
real, dimension (numst) :: masout_ri real, dimension (numst) :: masout_ri
...@@ -48,222 +66,438 @@ ...@@ -48,222 +66,438 @@
real, dimension (numst) :: masout_ct real, dimension (numst) :: masout_ct
real, dimension (numst) :: masout_ckt real, dimension (numst) :: masout_ckt
type (data_in) in type (meteoDataType) meteo_cell
type (data_outdef) out type (sfxDataType) sfx_cell
type (data_par) par1
type (data_lutyp) lu1 type (data_lutyp) lu1
integer i,numst integer i,numst
do i=1,numst do i = 1, numst
in%ws=mas1_w(i) meteo_cell = meteoDataType(&
in%dt=mas1_dt(i) h = meteo%h(i), &
in%st=mas1_st(i) U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
in%dq=mas1_dq(i) z0_m = meteo%z0_m(i))
in%cflh=mas1_cflh(i)
in%z0in=mas1_z0in(i) call surf_flux(meteo_cell, sfx_cell, numerics, lu1)
!if ((in%ws == in%ws).and.(in%dt == in%dt).and.(in%st == in%st).and.(in%dq == in%dq)) then masout_zl(i)=sfx_cell%zeta
CALL surf_flux(in, out, par1, lu1) masout_ri(i)=sfx_cell%Rib
!end if masout_re(i)=sfx_cell%Re
masout_lnzuzt(i)=sfx_cell%lnzuzt
masout_zl(i)=out%zl masout_zu(i)=sfx_cell%z0_m
masout_ri(i)=out%ri masout_ztout(i)=sfx_cell%z0_t
masout_re(i)=out%re masout_rith(i)=sfx_cell%Rib_conv_lim
masout_lnzuzt(i)=out%lnzuzt masout_cm(i)=sfx_cell%cm
masout_zu(i)=out%zu masout_ch(i)=sfx_cell%ch
masout_ztout(i)=out%ztout masout_ct(i)=sfx_cell%ct
masout_rith(i)=out%rith masout_ckt(i)=sfx_cell%ckt
masout_cm(i)=out%cm
masout_ch(i)=out%ch
masout_ct(i)=out%ct
masout_ckt(i)=out%ckt
enddo enddo
END SUBROUTINE surf_fluxMAS END SUBROUTINE surf_flux_vec
function f_m_conv(zeta)
real f_m_conv
real zeta
f_m_conv = (1.0 - alpha_m * zeta)**0.25
end function f_m_conv
function f_h_conv(zeta)
real f_h_conv
real zeta
f_h_conv = (1.0 - alpha_h * zeta)**0.5
end function f_h_conv
! stability functions (neutral)
! --------------------------------------------------------------------------------
subroutine get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
! ----------------------------------------------------------------------------
zeta = 0.0
psi_m = log(h0_m)
psi_h = log(h0_t) / Pr_t_0_inv
if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv
end subroutine
! --------------------------------------------------------------------------------
! stability functions (stable)
! --------------------------------------------------------------------------------
subroutine get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: Rib
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
real :: Rib_coeff
real :: psi0_m, psi0_h
real :: phi
real :: c
! ----------------------------------------------------------------------------
psi0_m = alog(h0_m)
psi0_h = B / psi0_m
Rib_coeff = beta_m * Rib
c = (psi0_h + 1.0) / Pr_t_0_inv - 2.0 * Rib_coeff
zeta = psi0_m * (sqrt(c**2 + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / &
(2.0 * beta_m * (1.0 - Rib_coeff))
phi = beta_m * zeta
psi_m = psi0_m + phi
psi_h = (psi0_m + B) / Pr_t_0_inv + phi
end subroutine
! --------------------------------------------------------------------------------
! stability functions (convection-semistrong)
! --------------------------------------------------------------------------------
subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: Rib
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
integer, intent(in) :: maxiters
real :: zeta0_m, zeta0_h
real :: x0, x1, y0, y1
integer :: i
! ----------------------------------------------------------------------------
psi_m = log(h0_m)
psi_h = log(h0_t)
if (abs(B) < 1.0e-10) psi_h = psi_m
zeta = Rib * Pr_t_0_inv * psi_m**2 / psi_h
do i = 1, maxiters + 1
zeta0_m = zeta / h0_m
zeta0_h = zeta / h0_t
if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
y1 = (1.0 - alpha_m * zeta)**0.25e0
x1 = sqrt(1.0 - alpha_h_fix * zeta)
y0 = (1.0 - alpha_m * zeta0_m)**0.25e0
x0 = sqrt(1.0 - alpha_h_fix * zeta0_h)
y0 = max(y0, 1.000001e0)
x0 = max(x0, 1.000001e0)
psi_m = log((y1 - 1.0e0)*(y0 + 1.0e0)/((y1 + 1.0e0)*(y0 - 1.0e0))) + 2.0e0*(atan(y1) - atan(y0))
psi_h = log((x1 - 1.0e0)*(x0 + 1.0e0)/((x1 + 1.0e0)*(x0 - 1.0e0))) / Pr_t_0_inv
if (i == maxiters + 1) exit
zeta = Rib * psi_m**2 / psi_h
end do
end subroutine
! --------------------------------------------------------------------------------
! stability functions (convection)
! --------------------------------------------------------------------------------
subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, &
zeta_conv_lim, x10, y10, &
h0_m, h0_t, B, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: Rib
real, intent(in) :: zeta_conv_lim
real, intent(in) :: x10, y10
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
integer, intent(in) :: maxiters
real :: zeta0_m, zeta0_h
real :: x0, y0
real :: p1, p0
real :: a1, b1, c, f
integer :: i
! ----------------------------------------------------------------------------
p1 = 2.0 * atan(y10) + log((y10 - 1.0) / (y10 + 1.0))
p0 = log((x10 - 1.0) / (x10 + 1.0))
zeta = zeta_conv_lim
do i = 1, maxiters + 1
zeta0_m = zeta / h0_m
zeta0_h = zeta / h0_t
if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
a1 = (zeta_conv_lim / zeta)**(1.0 / 3.0)
x0 = sqrt(1.0 - alpha_h_fix * zeta0_h)
y0 = (1.0 - alpha_m * zeta0_m)**0.25
c = log((x0 + 1.0)/(x0 - 1.0))
b1 = -2.0*atan(y0) + log((y0 + 1.0)/(y0 - 1.0))
f = 3.0 * (1.0 - a1)
psi_m = f / y10 + p1 + b1
psi_h = (f / x10 + p0 + c) / Pr_t_0_inv
if (i == maxiters + 1) exit
zeta = Rib * psi_m**2 / psi_h
end do
end subroutine
! --------------------------------------------------------------------------------
! define convection limit
! --------------------------------------------------------------------------------
subroutine get_convection_lim(zeta_lim, Rib_lim, x10, y10, &
h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: zeta_lim, Rib_lim
real, intent(out) :: x10, y10
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
real :: an4, an5
real :: psi_m, psi_h
! ----------------------------------------------------------------------------
an5=(Pr_t_inf_inv / Pr_t_0_inv)**4
zeta_lim = (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 = f_m_conv(zeta_lim)
x10 = f_h_conv(zeta_lim)
! ......definition of r-prim......
an4 = zeta_lim / h0_m
an5 = zeta_lim / h0_t
if (abs(B) < 1.0e-10) an5=an4
an5 = sqrt(1.0 - alpha_h_fix * an5)
an4 = (1.0 - alpha_m * an4)**0.25
psi_m = 2.0 * (atan(y10) - atan(an4)) + alog((y10 - 1.0) * (an4 + 1.0)/((y10 + 1.0) * (an4 - 1.0)))
psi_h = alog((x10 - 1.0) * (an5 + 1.0)/((x10 + 1.0) * (an5 - 1.0))) / Pr_t_0_inv
Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)
end subroutine
! --------------------------------------------------------------------------------
! charnock roughness
! --------------------------------------------------------------------------------
subroutine get_charnock_roughness(z0_m, u_dyn0, U, h, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_m, u_dyn0
real, intent(in) :: U, h
integer, intent(in) :: maxiters
real :: U10m
real :: a1, c1, y1, c1min, f
integer :: i, j
! ----------------------------------------------------------------------------
U10m = U
a1 = 0.0e0
y1 = 25.0e0
c1min = log(h_charnock) / kappa
do i = 1, maxiters
f = c1_charnock - 2.0 * log(U10m)
do j = 1, maxiters
c1 = (f + 2.0 * log(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 = log(1.0 + c2_charnock * ((y1 / U10m)**3)) / kappa
c1 = max(c1 - a1, c1min)
y1 = c1
end do
z0_m = h_charnock * exp(-c1 * kappa)
z0_m = max(z0_m,0.000015e0)
U10m = U*alog(h_charnock/z0_m)/(alog(h/z0_m))
end do
! --- define dynamic velocity in neutral conditions
u_dyn0 = U10m / c1
end subroutine
! --------------------------------------------------------------------------------
SUBROUTINE surf_flux(meteo, sfx, numerics, lu)
type (sfxDataType), intent(out) :: sfx
type (meteoDataType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
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_lutyp) lu
real ws, dt, dq, cflh, z0in real h ! surface flux layer height [m]
integer it
real z0_m, z0_t ! aerodynamic & thermal roughness [m]
real B ! = ln(z0_m / z0_t) [n/d]
real h0_m, h0_t ! = h / z0_m, h / z0_h [n/d]
real Re ! roughness Reynolds number = u_dyn0 * z0 / nu [n/d]
real u_dyn0 ! dynamic velocity in neutral conditions [m/s]
real zeta ! = z/L [n/d]
real zeta_conv_lim ! z/L critical value for matching free convection limit [n/d]
real Rib ! bulk Richardson number
real Rib_conv_lim ! Ri-bulk critical value for matching free convection limit [n/d]
real psi_m, psi_h ! universal functions [momentum] & [thermal]
real U ! wind velocity at h [m/s]
real Tsemi !
integer lu_indx 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 ! output ... !
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 an4, o, am
real d, c4, c1min, c0, c, b1, an0 real c4, c0, an0
! ...
! local unnamed !
real x10, y10
real al
real q4, t4
! .....
! just local variables
real f1, a1
! ....
integer is_ocean
integer i, j integer i, j
ws=in%ws
dt=in%dt Tsemi = meteo%Tsemi
Tsurf=in%st
dq=in%dq U = meteo%U
cflh=in%cflh t4 = meteo%dT
z0in=in%z0in q4 = meteo%dQ
it=par%it h = meteo%h
z0_m = meteo%z0_m
u=ws
t4=dt if(z0_m < 0.0) then
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 !> @brief Test - definition z0 of sea surface
!> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha) !> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
d0max = 8.0e0 is_ocean = 1
U10m=u
a1=0.0e0 call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
y1=25.0e0
c1min=alog(h1/1.0e0)/kappa ! --- define relative height
do i=1,it h0_m = h / z0_m
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 else
! ......parameters from viscosity sublayer...... ! ......parameters from viscosity sublayer......
h0=h/z0 is_ocean = 0
u3=u*kappa/alog(h0)
! --- define relative height
h0_m = h / z0_m
! --- define dynamic velocity in neutral conditions
u_dyn0 = U * kappa / log(h0_m)
end if
Re = u_dyn0 * z0_m / nu_air
if(Re <= Re_rough_min) then
B = B1_rough * alog(B3_rough * Re) + B2_rough
else
! B4 takes into account Re value at z' ~ O(10) z0
B = B4_rough * (Re**B2_rough)
end if end if
x7=u3*z0/an ! --- apply max restriction based on surface type
if(x7 <= x8) then if (is_ocean == 1) then
d0=an1*alog(al1*x7)+an2 B = min(B, B_max_ocean)
else else
d0=al2*(x7**0.45e0) B = min(B, B_max_land)
end if end if
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
! --- define relative height [thermal]
h0_t = h / z0_t
! ......humidity stratification and ri-number...... ! ......humidity stratification and ri-number......
al=g/Tsurf al = g / Tsemi
d0=min(d0,d0max) Rib = al * h * (t4 + 0.61e0 * Tsemi * q4) / U**2
Rib=al*h*(t4+0.61e0*Tsurf*q4)/u**2
d00=d0 ! --- define free convection transition zeta = z/L value
zt=z0/exp(d00) call get_convection_lim(zeta_conv_lim, Rib_conv_lim, x10, y10, &
h00=h/zt h0_m, h0_t, B)
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 if (Rib > 0.0e0) then
! ......stable stratification...... ! ......stable stratification......
!write (*,*) 'stable' !write (*,*) 'stable'
Rib=min(Rib,r0) Rib = min(Rib, Rib_max)
f=alog(h0) call get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
f1=d0/f
a1=b4*Rib f1 = beta_m * zeta
a2ch=(f1+1.0e0)/Pr_t_0_inv-2.0e0*a1 am = 1.0 + f1
d3=f*(sqrt(a2ch**2+4.0e0*a1*(1.0e0-a1))-a2ch)/(2.0e0*b4*(1.0e0-a1)) o = 1.0/Pr_t_0_inv + f1
f1=b4*d3
f4=f+f1 else if (Rib < Rib_conv_lim) then
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..... ! ......strong instability.....
!write (*,*) 'instability' !write (*,*) 'instability'
d3=d1 call get_psi_convection(psi_m, psi_h, zeta, Rib, &
do i=1,it+1 zeta_conv_lim, x10, y10, h0_m, h0_t, B, numerics%maxiters_convection)
d=d3/h0
dd=d3/h00 a1 = (zeta_conv_lim / zeta)**(1.0/3.0)
if (abs(d0) < 1.0e-10) dd=d am = a1 / y10
a1=(d1/d3)**(1.0e0/3.0e0) o = a1 / (Pr_t_0_inv * x10)
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 else if (Rib > -0.001e0) then
! ......nearly neutral...... ! ......nearly neutral......
write (*,*) 'neutral' write (*,*) 'neutral'
f4=alog(h0)
f0=ft0/Pr_t_0_inv call get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
if (abs(d0) < 1.0e-10) f0=f4/Pr_t_0_inv
am=1.0e0 am = 1.0e0
o=1.0e0/Pr_t_0_inv o = 1.0e0 / Pr_t_0_inv
else else
! ......week and semistrong instability...... ! ......week and semistrong instability......
!write (*,*) 'semistrong' !write (*,*) 'semistrong'
f1=alog(h0)
if (abs(d0) < 1.0e-10) ft0=f1 call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection)
d3=Rib*Pr_t_0_inv*f1**2/ft0
do i=1,it+1 am = (1.0 - alpha_m * zeta)**(-0.25)
d=d3/h0 o = 1.0/(Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta))
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 end if
! ......computation of cu,co,k(h),alft ! ......computation of cu,co,k(h),alft
c4=kappa/f4 c4=kappa/psi_m
c0=kappa/f0 c0=kappa/psi_h
an4=kappa*c4*u*h/am an4=kappa*c4*U*h/am
an0=am/o an0=am/o
! ......exit...... ! ......exit......
out%zl=d3 sfx%zeta = zeta
out%ri=Rib sfx%Rib = Rib
out%re=x7 sfx%Re = Re
out%lnzuzt=d00 sfx%lnzuzt=B
out%zu=z0 sfx%z0_m = z0_m
out%ztout=zt sfx%z0_t = z0_t
out%rith=r1 sfx%Rib_conv_lim = Rib_conv_lim
out%cm=c4 sfx%cm=c4
out%ch=c0 sfx%ch=c0
out%ct=an4 sfx%ct=an4
out%ckt=an0 sfx%ckt=an0
return return
END SUBROUTINE surf_flux END SUBROUTINE surf_flux
......
inputdata.f90
+ 6
− 6
View file @ c1e705f6
...@@ -2,12 +2,12 @@ module INPUTDATA ...@@ -2,12 +2,12 @@ module INPUTDATA
REAl, DIMENSION (6) :: AR1 REAl, DIMENSION (6) :: AR1
REAl, DIMENSION (11) :: AR2 REAl, DIMENSION (11) :: AR2
INTEGER, PARAMETER :: TEST=1 ! probably IT before renaming INTEGER, PARAMETER :: TEST=1 ! probably IT before renaming
INTEGER nums,ioer,mk !INTEGER nums,ioer,mk
REAL HFX, MFX, zL, betta !REAL HFX, MFX, zL, betta
REAL U, T4,c0,c4, T1,H,z0h !REAL U, T4,c0,c4, T1,H,z0h
REAL ws, deltaT, semisumT, D00, Z0, ZT, deltaQ !REAL ws, deltaT, semisumT, D00, Z0, ZT, deltaQ
REAL R6,R1 !REAL R6,R1
REAL AN5, D1, Y10, X10, P1, P0 !REAL AN5, Y10, X10, P1, P0
!C*==================================================================== !C*====================================================================
!C* .....DEFENITION OF DRAG AND HEAT EXCHANGE COEFFICIENTS...... = !C* .....DEFENITION OF DRAG AND HEAT EXCHANGE COEFFICIENTS...... =
......
main_drag.f90
+ 41
− 35
View file @ c1e705f6
...@@ -4,11 +4,11 @@ ...@@ -4,11 +4,11 @@
USE inputdata USE inputdata
USE drag3 USE drag3
type (data_in):: data_in1 type (meteoDataType):: data_in1
type (data_outdef) :: data_outdef1 type (sfxDataType) :: data_outdef1
type (data_par) :: data_par1 type (numericsType) :: data_par1
type (data_lutyp) :: data_lutyp1 type (data_lutyp) :: data_lutyp1
...@@ -18,15 +18,15 @@ ...@@ -18,15 +18,15 @@
character(len = 50) :: filename_out character(len = 50) :: filename_out
character(len = 50) :: filename_in2 character(len = 50) :: filename_in2
type :: datatype_inMAS1 !type :: datatype_inMAS1
real, allocatable :: mas_w(:) ! ! real, allocatable :: mas_w(:) !
real, allocatable :: mas_dt(:) ! real, allocatable :: mas_dt(:)
real, allocatable :: mas_st(:) ! real, allocatable :: mas_st(:)
real, allocatable :: mas_dq(:) ! real, allocatable :: mas_dq(:)
real, allocatable :: mas_cflh(:) ! real, allocatable :: mas_cflh(:)
real, allocatable :: mas_z0in(:) ! real, allocatable :: mas_z0in(:)
end type !end type
type(datatype_inMAS1) :: data_inMAS type(meteoDataVecType) :: meteo
!input !input
! mas_w - abs(wind velocity) at constant flux layer (cfl) hight (m/s) ! mas_w - abs(wind velocity) at constant flux layer (cfl) hight (m/s)
...@@ -86,19 +86,19 @@ ...@@ -86,19 +86,19 @@
open (2, file=filename_out) open (2, file=filename_out)
numst=0 numst=0
do WHILE (ioer.eq.0) do WHILE (ioer.eq.0)
read (1,*, iostat=ioer) data_in1%ws, data_in1%dt, data_in1%st, data_in1%dq read (1,*, iostat=ioer) data_in1%U, data_in1%dT, data_in1%Tsemi, data_in1%dQ
numst=numst+1 numst=numst+1
enddo enddo
close (1) close (1)
numst=numst-1 numst=numst-1
allocate(data_inMAS%mas_w(numst)) allocate(meteo%h(numst))
allocate(data_inMAS%mas_dt(numst)) allocate(meteo%U(numst))
allocate(data_inMAS%mas_st(numst)) allocate(meteo%dT(numst))
allocate(data_inMAS%mas_dq(numst)) allocate(meteo%Tsemi(numst))
allocate(data_inMAS%mas_cflh(numst)) allocate(meteo%dQ(numst))
allocate(data_inMAS%mas_z0in(numst)) allocate(meteo%z0_m(numst))
...@@ -118,23 +118,29 @@ ...@@ -118,23 +118,29 @@
open (1, file= filename_in, status ='old') open (1, file= filename_in, status ='old')
read (11, *) cflh, z0in read (11, *) cflh, z0in
do i=1,numst do i=1,numst
read (1,*) data_in1%ws, data_in1%dt, data_in1%st, data_in1%dq read (1,*) data_in1%U, data_in1%dT, data_in1%Tsemi, data_in1%dQ
data_inMAS%mas_w(i)=data_inMAS%mas_w(i)+data_in1%ws
data_inMAS%mas_dt(i)=data_inMAS%mas_dt(i)+data_in1%dt meteo%h(i)=cflh
data_inMAS%mas_st(i)=data_inMAS%mas_st(i)+data_in1%st meteo%U(i) = meteo%U(i)+data_in1%U
data_inMAS%mas_dq(i)=data_inMAS%mas_dq(i)+data_in1%dq meteo%dT(i) = meteo%dT(i)+data_in1%dT
data_inMAS%mas_cflh(i)=cflh meteo%Tsemi(i) = meteo%Tsemi(i)+data_in1%Tsemi
data_inMAS%mas_z0in(i)=z0in meteo%dQ(i) = meteo%dQ(i)+data_in1%dQ
meteo%z0_m(i)=z0in
enddo enddo
CALL surf_flux_vec(meteo, &
CALL surf_fluxMAS(data_inMAS%mas_w, data_inMAS%mas_dt, data_inMAS%mas_st, data_inMAS%mas_dq,&
data_inMAS%mas_cflh, data_inMAS%mas_z0in,&
data_outMAS%masout_zl, data_outMAS%masout_ri, data_outMAS%masout_re, data_outMAS%masout_lnzuzt,& data_outMAS%masout_zl, data_outMAS%masout_ri, data_outMAS%masout_re, data_outMAS%masout_lnzuzt,&
data_outMAS%masout_zu,data_outMAS%masout_ztout,data_outMAS%masout_rith,data_outMAS%masout_cm,& data_outMAS%masout_zu,data_outMAS%masout_ztout,data_outMAS%masout_rith,data_outMAS%masout_cm,&
data_outMAS%masout_ch,data_outMAS%masout_ct,data_outMAS%masout_ckt,& data_outMAS%masout_ch,data_outMAS%masout_ct,data_outMAS%masout_ckt,&
data_par1, data_lutyp1,numst) data_par1, data_lutyp1,numst)
!CALL surf_fluxMAS(meteo%mas_w, meteo%mas_dt, meteo%mas_st, meteo%mas_dq,&
! meteo%mas_cflh, meteo%mas_z0in,&
! data_outMAS%masout_zl, data_outMAS%masout_ri, data_outMAS%masout_re, data_outMAS%masout_lnzuzt,&
! data_outMAS%masout_zu,data_outMAS%masout_ztout,data_outMAS%masout_rith,data_outMAS%masout_cm,&
! data_outMAS%masout_ch,data_outMAS%masout_ct,data_outMAS%masout_ckt,&
! data_par1, data_lutyp1,numst)
do i=1,numst do i=1,numst
write (2,20) data_outMAS%masout_zl(i), data_outMAS%masout_ri(i), data_outMAS%masout_re(i), data_outMAS%masout_lnzuzt(i),& write (2,20) data_outMAS%masout_zl(i), data_outMAS%masout_ri(i), data_outMAS%masout_re(i), data_outMAS%masout_lnzuzt(i),&
data_outMAS%masout_zu(i), data_outMAS%masout_ztout(i), data_outMAS%masout_rith(i), data_outMAS%masout_cm(i),& data_outMAS%masout_zu(i), data_outMAS%masout_ztout(i), data_outMAS%masout_rith(i), data_outMAS%masout_cm(i),&
...@@ -142,12 +148,12 @@ ...@@ -142,12 +148,12 @@
enddo enddo
deallocate(data_inMAS%mas_w) deallocate(meteo%h)
deallocate(data_inMAS%mas_dt) deallocate(meteo%U)
deallocate(data_inMAS%mas_st) deallocate(meteo%dT)
deallocate(data_inMAS%mas_dq) deallocate(meteo%Tsemi)
deallocate(data_inMAS%mas_cflh) deallocate(meteo%dQ)
deallocate(data_inMAS%mas_z0in) deallocate(meteo%z0_m)
deallocate(data_outMAS%masout_zl) deallocate(data_outMAS%masout_zl)
......
param.f90
+ 39
− 25
View file @ c1e705f6
...@@ -2,40 +2,54 @@ module param ...@@ -2,40 +2,54 @@ module param
implicit none implicit none
! acceleration due to gravity [m/s^2] ! acceleration due to gravity [m/s^2]
real, parameter :: g = 9.81 real, parameter :: g = 9.81
! molecular Prandtl number (air) ! molecular Prandtl number (air) [n/d]
real, parameter :: Pr_m = 0.71 real, parameter :: Pr_m = 0.71
! kinematic viscosity of air [m^2/s]
real, parameter :: nu_air = 0.000015e0
! von Karman constant [n/d] ! von Karman constant [n/d]
real, parameter :: kappa = 0.40 real, parameter :: kappa = 0.40
! inverse Prandtl (turbulent) number in neutral conditions ! inverse Prandtl (turbulent) number in neutral conditions [n/d]
real, parameter :: Pr_t_0_inv = 1.15 real, parameter :: Pr_t_0_inv = 1.15
! inverse Prandtl (turbulent) number in free convection ! inverse Prandtl (turbulent) number in free convection [n/d]
real, parameter :: Pr_t_inf_inv = 3.5 real, parameter :: Pr_t_inf_inv = 3.5
! stability function coeff. [= g4 & g10 in deprecated code] ! stability function coeff. (unstable) [= g4 & g10 in deprecated code]
real, parameter :: alpha_m = 16.0 real, parameter :: alpha_m = 16.0
real, parameter :: alpha_h = 16.0 real, parameter :: alpha_h = 16.0
real, parameter :: alpha_h_fix = 1.2
! stability function coeff. (stable)
real, parameter :: beta_m = 4.7
real, parameter :: beta_h = beta_m
! --- max Ri-bulk value in stable case ( < 1 / beta_m )
real, parameter :: Rib_max = 0.9 / beta_m
! Re fully roughness minimum value [n/d]
real, parameter :: Re_rough_min = 16.3
! roughness model coeff. [n/d]
! --- transitional mode
! B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2
real, parameter :: B1_rough = 5.0 / 6.0
real, parameter :: B2_rough = 0.45
real, parameter :: B3_rough = kappa * Pr_m
! --- fully rough mode (Re > Re_rough_min)
! B = B4 * Re^(B2)
real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)
real, parameter :: B_max_land = 2.0
real, parameter :: B_max_ocean = 8.0
! Charnock parameters
! z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g)
real, parameter :: gamma_c = 0.0144
real, parameter :: Re_visc_min = 0.111
real, parameter :: h_charnock = 10.0
real, parameter :: c1_charnock = log(h_charnock * (g / gamma_c))
real, parameter :: c2_charnock = Re_visc_min * nu_air * c1_charnock
real, parameter :: b4=4.7e0
real, parameter :: alfam=.0144e0
real, parameter :: betam=.111e0
real, parameter :: an=.000015e0
real, parameter :: h1=10.0e0
real, parameter :: x8=16.3e0
real, parameter :: an1=5.0e0/6.0e0
real, parameter :: an2=.45e0
real, parameter :: al1=kappa*Pr_m
real, parameter :: g0=1.2
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
!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 ! AN5=(A6/A0)**4
! D1=(2.0E0*G10-AN5*G4-SQRT((AN5*G4)**2+4.0E0*AN5*G10*(G10-G4)))/(2.0E0*G10**2) ! 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 ! Y10=(1.0E0-G4*D1)**.25E0
......
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