From c1e705f6332b4743c682a294751254b81d16c28a Mon Sep 17 00:00:00 2001
From: Evgeny Mortikov <evgeny.mortikov@gmail.com>
Date: Sun, 17 Dec 2023 04:58:11 +0300
Subject: [PATCH] major code update
---
drag3.f90 | 624 ++++++++++++++++++++++++++++++++++----------------
inputdata.f90 | 12 +-
main_drag.f90 | 76 +++---
param.f90 | 64 ++++--
4 files changed, 515 insertions(+), 261 deletions(-)
diff --git a/drag3.f90 b/drag3.f90
index 6b20aeb..c93706f 100644
--- a/drag3.f90
+++ b/drag3.f90
@@ -1,21 +1,44 @@
MODULE drag3
USE param
- USE inputdata
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
- real, public:: ws, dt, st, dq, cflh, z0in
+ type, public :: sfxDataType
+ 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
-
- type, public:: data_outdef
- real, public:: zl, ri, re, lnzuzt, zu, ztout, rith, cm, ch, ct, ckt
+
+ type, public :: numericsType
+ 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
-
- type, public:: data_par
- integer, public :: it=10
- end type
type, public:: data_lutyp
integer, public :: lu_indx=1
@@ -23,18 +46,13 @@
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_flux_vec(meteo, &
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)
+ numerics, 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
+ type (meteoDataVecType), intent(in) :: meteo
+ type (numericsType), intent(in) :: numerics
real, dimension (numst) :: masout_zl
real, dimension (numst) :: masout_ri
@@ -48,222 +66,438 @@
real, dimension (numst) :: masout_ct
real, dimension (numst) :: masout_ckt
- type (data_in) in
- type (data_outdef) out
- type (data_par) par1
+ type (meteoDataType) meteo_cell
+ type (sfxDataType) sfx_cell
+
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
+ do i = 1, numst
+ meteo_cell = meteoDataType(&
+ h = meteo%h(i), &
+ U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
+ z0_m = meteo%z0_m(i))
+
+ call surf_flux(meteo_cell, sfx_cell, numerics, lu1)
+
+ masout_zl(i)=sfx_cell%zeta
+ masout_ri(i)=sfx_cell%Rib
+ masout_re(i)=sfx_cell%Re
+ masout_lnzuzt(i)=sfx_cell%lnzuzt
+ masout_zu(i)=sfx_cell%z0_m
+ masout_ztout(i)=sfx_cell%z0_t
+ masout_rith(i)=sfx_cell%Rib_conv_lim
+ masout_cm(i)=sfx_cell%cm
+ masout_ch(i)=sfx_cell%ch
+ masout_ct(i)=sfx_cell%ct
+ masout_ckt(i)=sfx_cell%ckt
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
- real ws, dt, dq, cflh, z0in
- integer it
+ real h ! surface flux layer height [m]
+
+ 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
- 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
+
+ ! output ... !
+ real an4, o, am
+ 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
- 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
+
+ Tsemi = meteo%Tsemi
+
+ U = meteo%U
+ t4 = meteo%dT
+ q4 = meteo%dQ
+ h = meteo%h
+ z0_m = meteo%z0_m
+
+ if(z0_m < 0.0) 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
+ is_ocean = 1
+
+ call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
+
+ ! --- define relative height
+ h0_m = h / z0_m
else
! ......parameters from viscosity sublayer......
- h0=h/z0
- u3=u*kappa/alog(h0)
+ is_ocean = 0
+
+ ! --- 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
- x7=u3*z0/an
- if(x7 <= x8) then
- d0=an1*alog(al1*x7)+an2
+ ! --- apply max restriction based on surface type
+ if (is_ocean == 1) then
+ B = min(B, B_max_ocean)
else
- d0=al2*(x7**0.45e0)
+ B = min(B, B_max_land)
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......
- 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......
+ al = g / Tsemi
+ Rib = al * h * (t4 + 0.61e0 * Tsemi * q4) / U**2
+
+ ! --- define free convection transition zeta = z/L value
+ call get_convection_lim(zeta_conv_lim, Rib_conv_lim, x10, y10, &
+ h0_m, h0_t, B)
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
+ Rib = min(Rib, Rib_max)
+ call get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
+
+ f1 = beta_m * zeta
+ am = 1.0 + f1
+ o = 1.0/Pr_t_0_inv + f1
+
+ else if (Rib < Rib_conv_lim) 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)
+ call get_psi_convection(psi_m, psi_h, zeta, Rib, &
+ zeta_conv_lim, x10, y10, h0_m, h0_t, B, numerics%maxiters_convection)
+
+ a1 = (zeta_conv_lim / zeta)**(1.0/3.0)
+ 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
+
+ call get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
+
+ 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))
+
+ call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection)
+
+ am = (1.0 - alpha_m * zeta)**(-0.25)
+ o = 1.0/(Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta))
end if
! ......computation of cu,co,k(h),alft
- c4=kappa/f4
- c0=kappa/f0
- an4=kappa*c4*u*h/am
+ c4=kappa/psi_m
+ c0=kappa/psi_h
+ 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
+ sfx%zeta = zeta
+ sfx%Rib = Rib
+ sfx%Re = Re
+ sfx%lnzuzt=B
+ sfx%z0_m = z0_m
+ sfx%z0_t = z0_t
+ sfx%Rib_conv_lim = Rib_conv_lim
+ sfx%cm=c4
+ sfx%ch=c0
+ sfx%ct=an4
+ sfx%ckt=an0
return
END SUBROUTINE surf_flux
diff --git a/inputdata.f90 b/inputdata.f90
index 3543420..1e9a006 100644
--- a/inputdata.f90
+++ b/inputdata.f90
@@ -2,12 +2,12 @@ module INPUTDATA
REAl, DIMENSION (6) :: AR1
REAl, DIMENSION (11) :: AR2
INTEGER, PARAMETER :: TEST=1 ! probably IT before renaming
- INTEGER nums,ioer,mk
- REAL HFX, MFX, zL, betta
- REAL U, T4,c0,c4, T1,H,z0h
- REAL ws, deltaT, semisumT, D00, Z0, ZT, deltaQ
- REAL R6,R1
- REAL AN5, D1, Y10, X10, P1, P0
+ !INTEGER nums,ioer,mk
+ !REAL HFX, MFX, zL, betta
+ !REAL U, T4,c0,c4, T1,H,z0h
+ !REAL ws, deltaT, semisumT, D00, Z0, ZT, deltaQ
+ !REAL R6,R1
+ !REAL AN5, Y10, X10, P1, P0
!C*====================================================================
!C* .....DEFENITION OF DRAG AND HEAT EXCHANGE COEFFICIENTS...... =
diff --git a/main_drag.f90 b/main_drag.f90
index 0d18fcf..776bee4 100644
--- a/main_drag.f90
+++ b/main_drag.f90
@@ -4,11 +4,11 @@
USE inputdata
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
@@ -18,15 +18,15 @@
character(len = 50) :: filename_out
character(len = 50) :: filename_in2
- type :: datatype_inMAS1
- real, allocatable :: mas_w(:) !
- real, allocatable :: mas_dt(:)
- real, allocatable :: mas_st(:)
- real, allocatable :: mas_dq(:)
- real, allocatable :: mas_cflh(:)
- real, allocatable :: mas_z0in(:)
- end type
- type(datatype_inMAS1) :: data_inMAS
+ !type :: datatype_inMAS1
+ ! real, allocatable :: mas_w(:) !
+ ! real, allocatable :: mas_dt(:)
+ ! real, allocatable :: mas_st(:)
+ ! real, allocatable :: mas_dq(:)
+ ! real, allocatable :: mas_cflh(:)
+ ! real, allocatable :: mas_z0in(:)
+ !end type
+ type(meteoDataVecType) :: meteo
!input
! mas_w - abs(wind velocity) at constant flux layer (cfl) hight (m/s)
@@ -86,19 +86,19 @@
open (2, file=filename_out)
numst=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
enddo
close (1)
numst=numst-1
- allocate(data_inMAS%mas_w(numst))
- allocate(data_inMAS%mas_dt(numst))
- allocate(data_inMAS%mas_st(numst))
- allocate(data_inMAS%mas_dq(numst))
- allocate(data_inMAS%mas_cflh(numst))
- allocate(data_inMAS%mas_z0in(numst))
+ allocate(meteo%h(numst))
+ allocate(meteo%U(numst))
+ allocate(meteo%dT(numst))
+ allocate(meteo%Tsemi(numst))
+ allocate(meteo%dQ(numst))
+ allocate(meteo%z0_m(numst))
@@ -118,23 +118,29 @@
open (1, file= filename_in, status ='old')
read (11, *) cflh, z0in
do i=1,numst
- read (1,*) data_in1%ws, data_in1%dt, data_in1%st, 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
- data_inMAS%mas_st(i)=data_inMAS%mas_st(i)+data_in1%st
- data_inMAS%mas_dq(i)=data_inMAS%mas_dq(i)+data_in1%dq
- data_inMAS%mas_cflh(i)=cflh
- data_inMAS%mas_z0in(i)=z0in
+ read (1,*) data_in1%U, data_in1%dT, data_in1%Tsemi, data_in1%dQ
+
+ meteo%h(i)=cflh
+ meteo%U(i) = meteo%U(i)+data_in1%U
+ meteo%dT(i) = meteo%dT(i)+data_in1%dT
+ meteo%Tsemi(i) = meteo%Tsemi(i)+data_in1%Tsemi
+ meteo%dQ(i) = meteo%dQ(i)+data_in1%dQ
+ meteo%z0_m(i)=z0in
enddo
-
- 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,&
+ CALL surf_flux_vec(meteo, &
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)
+ !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
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),&
@@ -142,12 +148,12 @@
enddo
- deallocate(data_inMAS%mas_w)
- deallocate(data_inMAS%mas_dt)
- deallocate(data_inMAS%mas_st)
- deallocate(data_inMAS%mas_dq)
- deallocate(data_inMAS%mas_cflh)
- deallocate(data_inMAS%mas_z0in)
+ deallocate(meteo%h)
+ deallocate(meteo%U)
+ deallocate(meteo%dT)
+ deallocate(meteo%Tsemi)
+ deallocate(meteo%dQ)
+ deallocate(meteo%z0_m)
deallocate(data_outMAS%masout_zl)
diff --git a/param.f90 b/param.f90
index 16cc630..4a1519e 100644
--- a/param.f90
+++ b/param.f90
@@ -2,40 +2,54 @@ module param
implicit none
! acceleration due to gravity [m/s^2]
real, parameter :: g = 9.81
- ! molecular Prandtl number (air)
+ ! molecular Prandtl number (air) [n/d]
real, parameter :: Pr_m = 0.71
+ ! kinematic viscosity of air [m^2/s]
+ real, parameter :: nu_air = 0.000015e0
! von Karman constant [n/d]
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
- ! inverse Prandtl (turbulent) number in free convection
+ ! inverse Prandtl (turbulent) number in free convection [n/d]
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_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
! 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
--
GitLab