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

major code update

parent c1e705f6
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
COMP_DRAG.txt
+ 1
− 1
View file @ 93ffc2f0
ifort -c inputdata.f90 fort -c inputdata.f90
ifort -c param.f90 ifort -c param.f90
ifort -c prmt.f90 ifort -c prmt.f90
ifort -c drag3.F ifort -c drag3.F
......
compile2.sh
+ 6
− 5
View file @ 93ffc2f0
#!/bin/bash #!/bin/bash
rm drag_ddt.exe *.o rm drag_ddt.exe *.o
gfortran -c -Wuninitialized inputdata.f90 gfortran -c -cpp -Wuninitialized inputdata.f90
gfortran -c -Wuninitialized param.f90 gfortran -c -cpp -Wuninitialized sfx_phys_const.f90
gfortran -c -Wuninitialized drag3.f90 gfortran -c -cpp -Wuninitialized sfx_esm_param.f90
gfortran -c -Wuninitialized main_drag.f90 gfortran -c -cpp -Wuninitialized sfx_esm.f90
gfortran -o drag_ddt.exe main_drag.o drag3.o inputdata.o param.o gfortran -c -Wuninitialized sfx_main.f90
gfortran -o drag_ddt.exe sfx_main.o sfx_esm.o inputdata.o sfx_esm_param.o sfx_phys_const.o
drag3.f90 deleted 100644 → 0
+ 0
− 505
View file @ c1e705f6
MODULE drag3
USE param
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 :: 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 :: 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_lutyp
integer, public :: lu_indx=1
end type
contains
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,&
numerics, lu1,numst)
type (meteoDataVecType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
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 (meteoDataType) meteo_cell
type (sfxDataType) sfx_cell
type (data_lutyp) lu1
integer i,numst
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_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
type (data_lutyp) lu
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
! 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
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)
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......
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
! --- apply max restriction based on surface type
if (is_ocean == 1) then
B = min(B, B_max_ocean)
else
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 / 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, 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'
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'
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'
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/psi_m
c0=kappa/psi_h
an4=kappa*c4*U*h/am
an0=am/o
! ......exit......
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
END MODULE drag3
\ No newline at end of file
makefile
+ 2
− 2
View file @ 93ffc2f0
RUN = drag.exe RUN = sfx
COMPILER ?= gnu COMPILER ?= gnu
FC_KEYS ?= FC_KEYS ?=
...@@ -10,7 +10,7 @@ ifeq ($(COMPILER),gnu) ...@@ -10,7 +10,7 @@ ifeq ($(COMPILER),gnu)
FC = gfortran FC = gfortran
endif endif
OBJ_F90 = inputdata.o param.o prmt.o OBJ_F90 = sfx_phys_const.o sfx_esm_param.o sfx_esm.o
OBJ_F = DRAG.o drag3.o OBJ_F = DRAG.o drag3.o
OBJ = $(OBJ_F90) $(OBJ_F) OBJ = $(OBJ_F90) $(OBJ_F)
......
sfx_esm.f90 0 → 100644
+ 593
− 0
View file @ 93ffc2f0
module sfx_esm
!> @brief main Earth System Model surface flux module
! macro defs.
! --------------------------------------------------------------------------------
#define SFX_FORCE_DEPRECATED_CODE
! --------------------------------------------------------------------------------
! modules used
! --------------------------------------------------------------------------------
use sfx_esm_param
! --------------------------------------------------------------------------------
! directives list
! --------------------------------------------------------------------------------
implicit none
private
! --------------------------------------------------------------------------------
! public interface
! --------------------------------------------------------------------------------
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
type, public :: meteoDataType
!> @brief meteorological input for surface flux calculation
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
!> @brief meteorological input for surface flux calculation
!> &details using arrays as input
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 :: sfxDataType
!> @brief surface flux output data
real :: zeta !> = z/L [n/d]
real :: Rib !> bulk Richardson number [n/d]
real :: Re !> Reynolds number [n/d]
real :: B !> = log(z0_m / z0_h) [n/d]
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 !> transfer coefficient for momentum [n/d]
real :: Ct !> transfer coefficient for heat [n/d]
real :: Km !> eddy viscosity coeff. at h [m^2/s]
real :: Pr_t_inv !> inverse turbulent Prandtl number at h [n/d]
end type
type, public :: sfxDataVecType
!> @brief surface flux output data
!> &details using arrays as output
real, allocatable :: zeta(:) !> = z/L [n/d]
real, allocatable :: Rib(:) !> bulk Richardson number [n/d]
real, allocatable :: Re(:) !> Reynolds number [n/d]
real, allocatable :: B(:) !> = log(z0_m / z0_h) [n/d]
real, allocatable :: z0_m(:) !> aerodynamic roughness [m]
real, allocatable :: z0_t(:) !> thermal roughness [m]
real, allocatable :: Rib_conv_lim(:) !> Ri-bulk convection critical value [n/d]
real, allocatable :: Cm(:) !> transfer coefficient for momentum [n/d]
real, allocatable :: Ct(:) !> transfer coefficient for heat [n/d]
real, allocatable :: Km(:) !> eddy viscosity coeff. at h [m^2/s]
real, allocatable :: Pr_t_inv(:) !> inverse turbulent Prandtl number at h [n/d]
end type
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
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
! --------------------------------------------------------------------------------
contains
! --------------------------------------------------------------------------------
subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
!> @brief surface flux calculation for array data
!> @details contains C/C++ & CUDA interface
! ----------------------------------------------------------------------------
type (sfxDataVecType), intent(inout) :: sfx
type (meteoDataVecType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
integer, intent(in) :: n
! ----------------------------------------------------------------------------
! --- local variables
type (meteoDataType) meteo_cell
type (sfxDataType) sfx_cell
integer i
! ----------------------------------------------------------------------------
do i = 1, n
#ifdef USE_DEPRECATED_CODE
#else
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 get_surface_fluxes(sfx_cell, meteo_cell, numerics)
call push_sfx_data(sfx, sfx_cell, i)
#endif
end do
end subroutine get_surface_fluxes_vec
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
subroutine push_sfx_data(sfx, sfx_cell, idx)
!> @brief helper subroutine for copying data in sfxDataVecType
! ----------------------------------------------------------------------------
type (sfxDataVecType), intent(inout) :: sfx
type (sfxDataType), intent(in) :: sfx_cell
integer, intent(in) :: idx
! ----------------------------------------------------------------------------
sfx%zeta(idx) = sfx_cell%zeta
sfx%Rib(idx) = sfx_cell%Rib
sfx%Re(idx) = sfx_cell%Re
sfx%B(idx) = sfx_cell%B
sfx%z0_m(idx) = sfx_cell%z0_m
sfx%z0_t(idx) = sfx_cell%z0_t
sfx%Rib_conv_lim(idx) = sfx_cell%Rib_conv_lim
sfx%Cm(idx) = sfx_cell%Cm
sfx%Ct(idx) = sfx_cell%Ct
sfx%Km(idx) = sfx_cell%Km
sfx%Pr_t_inv(idx) = sfx_cell%Pr_t_inv
end subroutine push_sfx_data
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
subroutine get_surface_fluxes(sfx, meteo, numerics)
!> @brief surface flux calculation for single cell
!> @details contains C/C++ interface
! ----------------------------------------------------------------------------
type (sfxDataType), intent(out) :: sfx
type (meteoDataType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
! ----------------------------------------------------------------------------
! --- meteo derived datatype name shadowing
! ----------------------------------------------------------------------------
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)
! ----------------------------------------------------------------------------
! --- local variables
! ----------------------------------------------------------------------------
real z0_t !> 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 u_dyn0 !> dynamic velocity in neutral conditions [m/s]
real Re !> roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
real zeta !> = z/L [n/d]
real Rib !> bulk Richardson number
real zeta_conv_lim !> z/L critical value for matching free convection limit [n/d]
real Rib_conv_lim !> Ri-bulk critical value for matching free convection limit [n/d]
real f_m_conv_lim !> stability function (momentum) value in free convection limit [n/d]
real f_h_conv_lim !> stability function (heat) value in free convection limit [n/d]
real psi_m, psi_h !> universal functions (momentum) & (heat) [n/d]
real phi_m, phi_h !> stability functions (momentum) & (heat) [n/d]
real Km !> eddy viscosity coeff. at h [m^2/s]
real Pr_t_inv !> invese Prandt number [n/d]
real Cm, Ct !> transfer coeff. for (momentum) & (heat) [n/d]
integer surface_type !> surface type = (ocean || land)
real fval !> just a shortcut for partial calculations
! ----------------------------------------------------------------------------
! --- shadowing names for clarity
U = meteo%U
Tsemi = meteo%Tsemi
dT = meteo%dT
dQ = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
! --- define surface type
if (z0_m < 0.0) then
surface_type = surface_ocean
else
surface_type = surface_land
end if
if (surface_type == surface_ocean) then
! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
! --- define relative height
h0_m = h / z0_m
endif
if (surface_type == surface_land) then
! --- define relative height
h0_m = h / z0_m
! --- define dynamic velocity in neutral conditions
u_dyn0 = U * kappa / log(h0_m)
end if
! --- define B = log(z0_m / z0_h)
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
! --- apply max restriction based on surface type
if (surface_type == surface_ocean) then
B = min(B, B_max_ocean)
endif
if (surface_type == surface_land) then
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
! --- define Ri-bulk
Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
! --- define free convection transition zeta = z/L value
call get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, &
h0_m, h0_t, B)
! --- get the fluxes
! ----------------------------------------------------------------------------
if (Rib > 0.0) then
! --- stable stratification block
! --- restrict bulk Ri value
! *: note that this value is written to output
Rib = min(Rib, Rib_max)
call get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
fval = beta_m * zeta
phi_m = 1.0 + fval
phi_h = 1.0/Pr_t_0_inv + fval
else if (Rib < Rib_conv_lim) then
! --- strong instability block
call get_psi_convection(psi_m, psi_h, zeta, Rib, &
zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, h0_m, h0_t, B, numerics%maxiters_convection)
fval = (zeta_conv_lim / zeta)**(1.0/3.0)
phi_m = fval / f_m_conv_lim
phi_h = fval / (Pr_t_0_inv * f_h_conv_lim)
else if (Rib > -0.001) then
! --- nearly neutral [-0.001, 0] block
call get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
phi_m = 1.0
phi_h = 1.0 / Pr_t_0_inv
else
! --- weak & semistrong instability block
call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection)
phi_m = (1.0 - alpha_m * zeta)**(-0.25)
phi_h = 1.0 / (Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta))
end if
! --- define transfer coeff. (momentum) & (heat)
Cm = kappa / psi_m
Ct = kappa / psi_h
! --- define eddy viscosity & inverse Prandtl number
Km = kappa * Cm * U * h / phi_m
Pr_t_inv = phi_m / phi_h
! --- setting output
sfx = sfxDataType(zeta = zeta, Rib = Rib, &
Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
Rib_conv_lim = Rib_conv_lim, &
Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
! convection universal functions shortcuts
! --------------------------------------------------------------------------------
function f_m_conv(zeta)
! ----------------------------------------------------------------------------
real :: f_m_conv
real, intent(in) :: 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, intent(in) :: zeta
! ----------------------------------------------------------------------------
f_h_conv = (1.0 - alpha_h * zeta)**0.5
end function f_h_conv
! --------------------------------------------------------------------------------
! universal functions
! --------------------------------------------------------------------------------
subroutine get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
!> @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !> universal functions
real, intent(out) :: zeta !> = z/L
real, intent(in) :: h0_m, h0_t !> = z/z0_m, z/z0_h
real, intent(in) :: B !> = log(z0_m / z0_h)
! ----------------------------------------------------------------------------
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
subroutine get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
!> @brief universal functions (momentum) & (heat): stable case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !> universal functions [n/d]
real, intent(out) :: zeta !> = z/L [n/d]
real, intent(in) :: Rib !> bulk Richardson number [n/d]
real, intent(in) :: h0_m, h0_t !> = z/z0_m, z/z0_h [n/d]
real, intent(in) :: B !> = log(z0_m / z0_h) [n/d]
! ----------------------------------------------------------------------------
! --- local variables
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
subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)
!> @brief universal functions (momentum) & (heat): semi-strong convection case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !> universal functions [n/d]
real, intent(out) :: zeta !> = z/L [n/d]
real, intent(in) :: Rib !> bulk Richardson number [n/d]
real, intent(in) :: h0_m, h0_t !> = z/z0_m, z/z0_h [n/d]
real, intent(in) :: B !> = log(z0_m / z0_h) [n/d]
integer, intent(in) :: maxiters !> maximum number of iterations
! --- local variables
real :: zeta0_m, zeta0_h
real :: f0_m, f0_h
real :: f_m, f_h
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
f_m = (1.0 - alpha_m * zeta)**0.25e0
f_h = sqrt(1.0 - alpha_h_fix * zeta)
f0_m = (1.0 - alpha_m * zeta0_m)**0.25e0
f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h)
f0_m = max(f0_m, 1.000001e0)
f0_h = max(f0_h, 1.000001e0)
psi_m = log((f_m - 1.0e0)*(f0_m + 1.0e0)/((f_m + 1.0e0)*(f0_m - 1.0e0))) + 2.0e0*(atan(f_m) - atan(f0_m))
psi_h = log((f_h - 1.0e0)*(f0_h + 1.0e0)/((f_h + 1.0e0)*(f0_h - 1.0e0))) / Pr_t_0_inv
! *: don't update zeta = z/L at last iteration
if (i == maxiters + 1) exit
zeta = Rib * psi_m**2 / psi_h
end do
end subroutine
subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, &
zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, &
h0_m, h0_t, B, maxiters)
!> @brief universal functions (momentum) & (heat): fully convective case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !> universal functions [n/d]
real, intent(out) :: zeta !> = z/L [n/d]
real, intent(in) :: Rib !> bulk Richardson number [n/d]
real, intent(in) :: h0_m, h0_t !> = z/z0_m, z/z0_h [n/d]
real, intent(in) :: B !> = log(z0_m / z0_h) [n/d]
integer, intent(in) :: maxiters !> maximum number of iterations
real, intent(in) :: zeta_conv_lim !> convective limit zeta
real, intent(in) :: f_m_conv_lim, f_h_conv_lim !> universal function shortcuts in limiting case
! ----------------------------------------------------------------------------
! --- local variables
real :: zeta0_m, zeta0_h
real :: f0_m, f0_h
real :: p_m, p_h
real :: a_m, a_h
real :: c_lim, f
integer :: i
! ----------------------------------------------------------------------------
p_m = 2.0 * atan(f_m_conv_lim) + log((f_m_conv_lim - 1.0) / (f_m_conv_lim + 1.0))
p_h = log((f_h_conv_lim - 1.0) / (f_h_conv_lim + 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
f0_m = (1.0 - alpha_m * zeta0_m)**0.25
f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h)
a_m = -2.0*atan(f0_m) + log((f0_m + 1.0)/(f0_m - 1.0))
a_h = log((f0_h + 1.0)/(f0_h - 1.0))
c_lim = (zeta_conv_lim / zeta)**(1.0 / 3.0)
f = 3.0 * (1.0 - c_lim)
psi_m = f / f_m_conv_lim + p_m + a_m
psi_h = (f / f_h_conv_lim + p_h + a_h) / Pr_t_0_inv
! *: don't update zeta = z/L at last iteration
if (i == maxiters + 1) exit
zeta = Rib * psi_m**2 / psi_h
end do
end subroutine
! --------------------------------------------------------------------------------
! convection limit definition
! --------------------------------------------------------------------------------
subroutine get_convection_lim(zeta_lim, Rib_lim, f_m_lim, f_h_lim, &
h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: zeta_lim !> limiting value of z/L
real, intent(out) :: Rib_lim !> limiting value of Ri-bulk
real, intent(out) :: f_m_lim, f_h_lim !> limiting values of universal functions shortcuts
real, intent(in) :: h0_m, h0_t !> = z/z0_m, z/z0_h [n/d]
real, intent(in) :: B !> = log(z0_m / z0_h) [n/d]
! ----------------------------------------------------------------------------
! --- local variables
real :: psi_m, psi_h
real :: f_m, f_h
real :: c
! ----------------------------------------------------------------------------
! --- define limiting value of zeta = z / L
c = (Pr_t_inf_inv / Pr_t_0_inv)**4
zeta_lim = (2.0 * alpha_h - c * alpha_m - &
sqrt((c * alpha_m)**2 + 4.0 * c * alpha_h * (alpha_h - alpha_m))) / (2.0 * alpha_h**2)
f_m_lim = f_m_conv(zeta_lim)
f_h_lim = f_h_conv(zeta_lim)
! --- universal functions
f_m = zeta_lim / h0_m
f_h = zeta_lim / h0_t
if (abs(B) < 1.0e-10) f_h = f_m
f_m = (1.0 - alpha_m * f_m)**0.25
f_h = sqrt(1.0 - alpha_h_fix * f_h)
psi_m = 2.0 * (atan(f_m_lim) - atan(f_m)) + alog((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0)))
psi_h = alog((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / Pr_t_0_inv
! --- bulk Richardson number
Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)
end subroutine
! --------------------------------------------------------------------------------
! charnock roughness definition
! --------------------------------------------------------------------------------
subroutine get_charnock_roughness(z0_m, u_dyn0, U, h, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_m !> aerodynamic roughness [m]
real, intent(out) :: u_dyn0 !> dynamic velocity in neutral conditions [m/s]
real, intent(in) :: h !> constant flux layer height [m]
real, intent(in) :: U !> abs(wind speed) [m/s]
integer, intent(in) :: maxiters !> maximum number of iterations
! ----------------------------------------------------------------------------
! --- local variables
real :: Uc ! wind speed at h_charnock [m/s]
real :: a, b, c, c_min
real :: f
integer :: i, j
! ----------------------------------------------------------------------------
Uc = U
a = 0.0
b = 25.0
c_min = log(h_charnock) / kappa
do i = 1, maxiters
f = c1_charnock - 2.0 * log(Uc)
do j = 1, maxiters
c = (f + 2.0 * log(b)) / 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) a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
c = max(c - a, c_min)
b = c
end do
z0_m = h_charnock * exp(-c * kappa)
z0_m = max(z0_m, 0.000015e0)
Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
end do
! --- define dynamic velocity in neutral conditions
u_dyn0 = Uc / c
end subroutine
! --------------------------------------------------------------------------------
end module sfx_esm
\ No newline at end of file
param.f90 → sfx_esm_param.f90
+ 71
− 0
View file @ 93ffc2f0
module param module sfx_esm_param
!> @brief ESM surface flux model parameters
!> @details all in SI units
! modules used
! --------------------------------------------------------------------------------
use sfx_phys_const
! --------------------------------------------------------------------------------
! directives list
! --------------------------------------------------------------------------------
implicit none implicit none
! acceleration due to gravity [m/s^2] ! --------------------------------------------------------------------------------
real, parameter :: g = 9.81
! molecular Prandtl number (air) [n/d] integer, public, parameter :: surface_ocean = 0 !> ocean surface
real, parameter :: Pr_m = 0.71 integer, public, parameter :: surface_land = 1 !> land surface
! 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 [n/d] !> 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 [n/d] !> 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. (unstable) [= 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 real, parameter :: alpha_h_fix = 1.2
! stability function coeff. (stable) !> stability function coeff. (stable)
real, parameter :: beta_m = 4.7 real, parameter :: beta_m = 4.7
real, parameter :: beta_h = beta_m real, parameter :: beta_h = beta_m
! --- max Ri-bulk value in stable case ( < 1 / beta_m ) !> --- max Ri-bulk value in stable case ( < 1 / beta_m )
real, parameter :: Rib_max = 0.9 / beta_m real, parameter :: Rib_max = 0.9 / beta_m
! Re fully roughness minimum value [n/d] !> Re fully roughness minimum value [n/d]
real, parameter :: Re_rough_min = 16.3 real, parameter :: Re_rough_min = 16.3
! roughness model coeff. [n/d] !> roughness model coeff. [n/d]
! --- transitional mode !> --- transitional mode
! B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2 !> B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2
real, parameter :: B1_rough = 5.0 / 6.0 real, parameter :: B1_rough = 5.0 / 6.0
real, parameter :: B2_rough = 0.45 real, parameter :: B2_rough = 0.45
real, parameter :: B3_rough = kappa * Pr_m real, parameter :: B3_rough = kappa * Pr_m
! --- fully rough mode (Re > Re_rough_min) !> --- fully rough mode (Re > Re_rough_min)
! B = B4 * Re^(B2) !> B = B4 * Re^(B2)
real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8) 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_land = 2.0
real, parameter :: B_max_ocean = 8.0 real, parameter :: B_max_ocean = 8.0
! Charnock parameters !> Charnock parameters
! z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g) !> z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g)
real, parameter :: gamma_c = 0.0144 real, parameter :: gamma_c = 0.0144
real, parameter :: Re_visc_min = 0.111 real, parameter :: Re_visc_min = 0.111
...@@ -50,14 +59,7 @@ module param ...@@ -50,14 +59,7 @@ module param
real, parameter :: c1_charnock = log(h_charnock * (g / gamma_c)) real, parameter :: c1_charnock = log(h_charnock * (g / gamma_c))
real, parameter :: c2_charnock = Re_visc_min * nu_air * c1_charnock real, parameter :: c2_charnock = Re_visc_min * nu_air * c1_charnock
! AN5=(A6/A0)**4 end module sfx_esm_param
! 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))
end module PARAM
......
main_drag.f90 → sfx_main.f90
+ 49
− 51
View file @ 93ffc2f0
PROGRAM main_ddt PROGRAM main_ddt
USE param use sfx_phys_const
USE sfx_esm_param
USE inputdata USE inputdata
USE drag3 USE sfx_esm
type (meteoDataType):: data_in1 type (meteoDataType):: data_in1
...@@ -11,7 +12,6 @@ ...@@ -11,7 +12,6 @@
type (numericsType) :: data_par1 type (numericsType) :: data_par1
type (data_lutyp) :: data_lutyp1
integer :: numst, i integer :: numst, i
real :: cflh, z0in real :: cflh, z0in
character(len = 50) :: filename_in character(len = 50) :: filename_in
...@@ -39,21 +39,21 @@ ...@@ -39,21 +39,21 @@
! lu_indx - 1 for land, 2 for sea, 3 for lake ! lu_indx - 1 for land, 2 for sea, 3 for lake
! test - file input ! test - file input
type :: datatype_outMAS1 !type :: datatype_outMAS1
real, allocatable :: masout_zl(:) ! real, allocatable :: masout_zl(:)
real, allocatable :: masout_ri(:) ! real, allocatable :: masout_ri(:)
real, allocatable :: masout_re(:) ! real, allocatable :: masout_re(:)
real, allocatable :: masout_lnzuzt(:) ! real, allocatable :: masout_lnzuzt(:)
real, allocatable :: masout_zu(:) ! real, allocatable :: masout_zu(:)
real, allocatable :: masout_ztout(:) ! real, allocatable :: masout_ztout(:)
real, allocatable :: masout_rith(:) ! real, allocatable :: masout_rith(:)
real, allocatable :: masout_cm(:) ! real, allocatable :: masout_cm(:)
real, allocatable :: masout_ch(:) ! real, allocatable :: masout_ch(:)
real, allocatable :: masout_ct(:) ! real, allocatable :: masout_ct(:)
real, allocatable :: masout_ckt(:) ! real, allocatable :: masout_ckt(:)
end type !end type
type(datatype_outMAS1) :: data_outMAS type(sfxDataVecType) :: data_outMAS
!output !output
!masout_zl - non-dimensional cfl hight !masout_zl - non-dimensional cfl hight
...@@ -101,18 +101,17 @@ ...@@ -101,18 +101,17 @@
allocate(meteo%z0_m(numst)) allocate(meteo%z0_m(numst))
allocate(data_outMAS%zeta(numst))
allocate(data_outMAS%masout_zl(numst)) allocate(data_outMAS%Rib(numst))
allocate(data_outMAS%masout_ri(numst)) allocate(data_outMAS%Re(numst))
allocate(data_outMAS%masout_re(numst)) allocate(data_outMAS%B(numst))
allocate(data_outMAS%masout_lnzuzt(numst)) allocate(data_outMAS%z0_m(numst))
allocate(data_outMAS%masout_zu(numst)) allocate(data_outMAS%z0_t(numst))
allocate(data_outMAS%masout_ztout(numst)) allocate(data_outMAS%Rib_conv_lim(numst))
allocate(data_outMAS%masout_rith(numst)) allocate(data_outMAS%Cm(numst))
allocate(data_outMAS%masout_cm(numst)) allocate(data_outMAS%Ct(numst))
allocate(data_outMAS%masout_ch(numst)) allocate(data_outMAS%Km(numst))
allocate(data_outMAS%masout_ct(numst)) allocate(data_outMAS%Pr_t_inv(numst))
allocate(data_outMAS%masout_ckt(numst))
open (11, file=filename_in2, status ='old') open (11, file=filename_in2, status ='old')
open (1, file= filename_in, status ='old') open (1, file= filename_in, status ='old')
...@@ -128,11 +127,11 @@ ...@@ -128,11 +127,11 @@
meteo%z0_m(i)=z0in meteo%z0_m(i)=z0in
enddo enddo
CALL surf_flux_vec(meteo, & CALL get_surface_fluxes_vec(data_outMAS, meteo, &
data_outMAS%masout_zl, data_outMAS%masout_ri, data_outMAS%masout_re, data_outMAS%masout_lnzuzt,& !data_outMAS%zeta, data_outMAS%Rib, data_outMAS%Re, data_outMAS%B,&
data_outMAS%masout_zu,data_outMAS%masout_ztout,data_outMAS%masout_rith,data_outMAS%masout_cm,& !data_outMAS%z0_m,data_outMAS%z0_t,data_outMAS%Rib_conv_lim,data_outMAS%Cm,&
data_outMAS%masout_ch,data_outMAS%masout_ct,data_outMAS%masout_ckt,& !data_outMAS%Ct,data_outMAS%Km,data_outMAS%Pr_t_inv,&
data_par1, data_lutyp1,numst) data_par1, numst)
!CALL surf_fluxMAS(meteo%mas_w, meteo%mas_dt, meteo%mas_st, meteo%mas_dq,& !CALL surf_fluxMAS(meteo%mas_w, meteo%mas_dt, meteo%mas_st, meteo%mas_dq,&
! meteo%mas_cflh, meteo%mas_z0in,& ! meteo%mas_cflh, meteo%mas_z0in,&
...@@ -142,9 +141,9 @@ ...@@ -142,9 +141,9 @@
! data_par1, data_lutyp1,numst) ! 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%zeta(i), data_outMAS%Rib(i), data_outMAS%Re(i), data_outMAS%B(i),&
data_outMAS%masout_zu(i), data_outMAS%masout_ztout(i), data_outMAS%masout_rith(i), data_outMAS%masout_cm(i),& data_outMAS%z0_m(i), data_outMAS%z0_t(i), data_outMAS%Rib_conv_lim(i), data_outMAS%Cm(i),&
data_outMAS%masout_ch(i), data_outMAS%masout_ct(i), data_outMAS%masout_ckt(i) data_outMAS%Ct(i), data_outMAS%Km(i), data_outMAS%Pr_t_inv(i)
enddo enddo
...@@ -155,18 +154,17 @@ ...@@ -155,18 +154,17 @@
deallocate(meteo%dQ) deallocate(meteo%dQ)
deallocate(meteo%z0_m) deallocate(meteo%z0_m)
deallocate(data_outMAS%zeta)
deallocate(data_outMAS%masout_zl) deallocate(data_outMAS%Rib)
deallocate(data_outMAS%masout_ri) deallocate(data_outMAS%Re)
deallocate(data_outMAS%masout_re) deallocate(data_outMAS%B)
deallocate(data_outMAS%masout_lnzuzt) deallocate(data_outMAS%z0_m)
deallocate(data_outMAS%masout_zu) deallocate(data_outMAS%z0_t)
deallocate(data_outMAS%masout_ztout) deallocate(data_outMAS%Rib_conv_lim)
deallocate(data_outMAS%masout_rith) deallocate(data_outMAS%Cm)
deallocate(data_outMAS%masout_cm) deallocate(data_outMAS%Ct)
deallocate(data_outMAS%masout_ch) deallocate(data_outMAS%Km)
deallocate(data_outMAS%masout_ct) deallocate(data_outMAS%Pr_t_inv)
deallocate(data_outMAS%masout_ckt)
10 format (f8.4,2x,f8.4) 10 format (f8.4,2x,f8.4)
......
sfx_phys_const.f90 0 → 100644
+ 12
− 0
View file @ 93ffc2f0
module sfx_phys_const
!> @brief general physics constants
!> @details all in SI units
implicit none
public
real, parameter :: g = 9.81 !> acceleration due to gravity [m/s^2]
real, parameter :: nu_air = 0.000015e0 !> kinematic viscosity of air [m^2/s]
real, parameter :: Pr_m = 0.71 !> molecular Prandtl number (air) [n/d]
end module sfx_phys_const
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment