Skip to content
Snippets Groups Projects
Commit 393b2460 authored by Anna Shestakova's avatar Anna Shestakova
Browse files

diag module for interpolation in surface layer

parent 77236674
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 7
− 2
View file @ 393b2460
...@@ -93,6 +93,11 @@ set(SOURCES_F ...@@ -93,6 +93,11 @@ set(SOURCES_F
srcF/sfx_api_inmcm.f90 srcF/sfx_api_inmcm.f90
srcF/sfx_api_term.f90 srcF/sfx_api_term.f90
) )
set(DIAG
diag/pbldia_new_sfx.f90
)
set(MAIN_F set(MAIN_F
srcF/sfx_main.f90) srcF/sfx_main.f90)
...@@ -162,11 +167,11 @@ endif(USE_CXX OR INCLUDE_CUDA) ...@@ -162,11 +167,11 @@ endif(USE_CXX OR INCLUDE_CUDA)
set(SOURCES set(SOURCES
${MEMPROC_HEADERS_CU} ${MEMPROC_SOURCES_CU} ${MEMPROC_HEADERS_CXX} ${MEMPROC_SOURCES_CXX} ${MEMPROC_HEADERS_CU} ${MEMPROC_SOURCES_CU} ${MEMPROC_HEADERS_CXX} ${MEMPROC_SOURCES_CXX}
${HEADERS_CU} ${SOURCES_CU} ${HEADERS_C} ${HEADERS_CXX} ${SOURCES_CXX} ${SOURCES_C} ${HEADERS_CU} ${SOURCES_CU} ${HEADERS_C} ${HEADERS_CXX} ${SOURCES_CXX} ${SOURCES_C}
${HEADERS_F} ${SOURCES_F} ${MAIN_F}) ${HEADERS_F} ${SOURCES_F} ${DIAG} ${MAIN_F})
set(SOURCES_LIB set(SOURCES_LIB
${MEMPROC_HEADERS_CU} ${MEMPROC_SOURCES_CU} ${MEMPROC_HEADERS_CXX} ${MEMPROC_SOURCES_CXX} ${MEMPROC_HEADERS_CU} ${MEMPROC_SOURCES_CU} ${MEMPROC_HEADERS_CXX} ${MEMPROC_SOURCES_CXX}
${HEADERS_CU} ${SOURCES_CU} ${HEADERS_CXX} ${HEADERS_CU} ${SOURCES_CU} ${HEADERS_CXX}
${SOURCES_CXX} ${SOURCES_C} ${HEADERS_F} ${SOURCES_F}) ${SOURCES_CXX} ${SOURCES_C} ${HEADERS_F} ${SOURCES_F} ${DIAG})
set(CMAKE_Fortran_FLAGS " -cpp ") set(CMAKE_Fortran_FLAGS " -cpp ")
......
diag/compile_simple.sh 0 → 100755
+ 21
− 0
View file @ 393b2460
#!/bin/bash
gfortran -c -cpp -Wuninitialized ../srcF/sfx_phys_const.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_esm_param.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_sheba_param.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_most_param.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_log_param.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_data.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_surface.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_esm.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_sheba.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_most.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_log.f90
gfortran -c -cpp -Wuninitialized pbldia_new_sfx.f90
gfortran -c -cpp -Wuninitialized diag_sfc_simple.f90
gfortran -o a.exe sfx_phys_const.o sfx_sheba_param.o sfx_esm_param.o sfx_most_param.o sfx_log_param.o sfx_data.o sfx_surface.o sfx_sheba.o sfx_esm.o sfx_most.o sfx_log.o pbldia_new_sfx.o diag_sfc_simple.o
rm *.o *.mod
diag/diag_sfc_simple.f90 0 → 100755
+ 192
− 0
View file @ 393b2460
! simple driver for wind, temperature and humidity interpolation in surface layer via pbldia
#include "../includeF/sfx_def.fi"
#include "sfx_diag_def.fi"
use sfx_data
use sfx_sheba, only: &
get_surface_fluxes_vec_sheba => get_surface_fluxes_vec, &
numericsType_sheba => numericsType
use sfx_esm, only: &
get_surface_fluxes_vec_esm => get_surface_fluxes_vec, &
numericsType_esm => numericsType
use sfx_most, only: &
get_surface_fluxes_vec_most => get_surface_fluxes_vec, &
numericsType_most => numericsType
use sfx_log, only: &
get_surface_fluxes_vec_log => get_surface_fluxes_vec, &
numericsType_log => numericsType
use pbldia_new_sfx
use sfx_phys_const
! input/output data
! --------------------------------------------------------------------------------
type(meteoDataVecType) :: meteo !< meteorological data (input)
type(meteoDataType) :: meteo_cell
type(sfxDataVecType) :: sfx !< surface fluxes (output)
#ifdef SHEBA
type(numericsType_sheba) :: numerics_sheba !< surface flux module (SHEBA) numerics parameters
#endif
#ifdef ESM
type(numericsType_esm) :: numerics_esm !< surface flux module (ESM) numerics parameters
#endif
#ifdef MOST
type(numericsType_most) :: numerics_most !< surface flux module (MOST) numerics parameters
#endif
#ifdef LOG
type(numericsType_log) :: numerics_log !< surface flux module (LOG) numerics parameters
#endif
integer :: num !< number of 'cells' in input
real,parameter::hmin=10.,hstep=10.,hmax=100. ! minimum, maximum heights and step (m) of desired wind profiles
real,parameter::undef=-999.0,R=287.,appa=0.287
real AR2(11),ARDIN(6),ARDOUT(3)
real z0m,z0h,cm,ct,dt,dq,h,zw,zt,L
real t1,u1,q1,th1,p1,Tin,Thin,u,Ts,Ths,Qin,Qs,Pin,Ps,z,Hpbl,Upbl,Tpbl,Qpbl
character(10) surf_type
logical pot_temp,calc_ps
!real wind_profile(int((hmax-hmin)/hstep)),temp_profile(int((hmax-hmin)/hstep)),q_profile(int((hmax-hmin)/hstep))
integer j,k
! user-defined input parameters for wind, temperature and humidity interpolation in the surface layer
!mandatory:
u = 5.0
Tin = 288.0
Ts = 283.0
Qin = 0.005
Qs = 0.008
z0m = 0.03 !momentum roughness length
pot_temp = .true. !if inout temperatures are potential
h = 25. ! input height
zw = 10 ! height to which wind interpolation is performed
zt = 2 ! height to wich temperature and humidity interpolation is performed
! auxiliary:
Pin = 1000. ! input pressure (hPa)
Ps = undef ! surface pressure (in hPa)
z0h = undef !thermal roughness length
cm = undef !transfer coefficient for momentum
ct = undef ! transfer coefficient for heat and humidity
L = undef ! MO length
surf_type='land' !needed for superadiabatic gradients correction
!----------end of user-defined variables
if(Ps == undef)then
calc_ps = .true.
endif
! check if temperature is in K
if(Tin.lt.200) Tin = Tin + 273.15
if(Ts.lt.200) Ts = Ts + 273.15
!check if input temperatures are potential or absolute
if(.not.pot_temp)then
Thin = Tin * (1000./Pin)**appa
if(calc_ps)then
Ps = Pin / exp(-g * h / (R * 0.5 * (Tin * (1 + 0.608 * Qin) + Ts * (1 + 0.608 * Qs))))
endif
Ths = Ts * (1000./Ps)**appa
else
Thin = Tin
Tin = Thin / ((1000./Pin)**appa)
if(calc_ps)then
Ps = Pin / exp(-g * h / (R * 0.5 * (Tin * (1 + 0.608 * Qin) + Ts * (1 + 0.608 * Qs))))
endif
Ths = Ts
Ts = Ths / ((1000./Ps)**appa)
endif
dt = Thin - Ths
dq = Qin - Qs
!calculation of Cm, Ct and L if they are not known
if(cm == undef.or.ct == undef.or.L == undef)then
!< @brief allocate input & output data
call allocate_meteo_vec(meteo, 1)
call allocate_sfx_vec(sfx, 1)
!< @brief read input data
meteo%h = h
meteo%U = u
meteo%dT = dt
meteo%Tsemi = 0.5*(Thin + Ths)
meteo%dQ = dq
meteo%z0_m = z0m
!< @brief calling flux module
#ifdef SHEBA
call get_surface_fluxes_vec_sheba(sfx, meteo, numerics_sheba, 1)
#endif
#ifdef ESM
call get_surface_fluxes_vec_esm(sfx, meteo, numerics_esm, 1)
#endif
#ifdef MOST
call get_surface_fluxes_vec_most(sfx, meteo, numerics_most, 1)
#endif
#ifdef LOG
call get_surface_fluxes_vec_log(sfx, meteo, numerics_log, 1)
#endif
ar2(1) = sfx%zeta(1)
ar2(5) = sfx%z0_m(1)
ar2(6) = sfx%z0_t(1)
ar2(8) = sfx%Cm(1)
ar2(9) = sfx%Ct(1)
else
ar2(1) = h / L
ar2(5) = z0m
ar2(6) = z0h
ar2(8) = cm
ar2(9) = ct
endif
ardin(1) = u
ardin(2) = dt
ardin(3) = dq
ardin(4) = h
ardin(5) = zw
ardin(6) = zt
!ardin(7) = Ths
!ardin(8) = Qs
!ardin(9) = Pin
!ardin(10) = Ps
!-----------------
#ifdef SHEBA
call pbldia_new_sheba(ar2(:),ardin(:),ardout(:))
#endif
#ifdef MOST
call pbldia_new_most(ar2(:),ardin(:),ardout(:))
#endif
#ifdef ESM
call pbldia_new_esm(ar2(:),ardin(:),ardout(:))
#endif
print *, "wind speed at z=",h," ",u,"; at z=",zw," ",ardout(1),"; at surface",0.
print *, "air temperature at z=",h," ",Tin,"; at z=",zt," ",ardout(2)+Ths,"; at surface",Ths
print *, "air humidity at z=",h," ",Qin,"; at z=",zt," ",ardout(3)+Qs,"; at surface",Qs
call deallocate_meteo_vec(meteo)
call deallocate_sfx_vec(sfx)
end
\ No newline at end of file
diag/pbldia_new_sfx.f90 0 → 100755
+ 252
− 0
View file @ 393b2460
module pbldia_new_sfx
use sfx_phys_const
implicit none
real,parameter,private::undef=-999.0,R=287.,appa=0.287
contains
! in this subroutine interpolation is performed using integration of universal functions from z1 to z2,
! where z1 is measurement or lowest model level height and z2 is required height
! The input-output structure is preserved as much as possible as in the climate model (pbldia subroutine in pblflx)
!C* DIAGNOSIS OF WIND VELOCITY (AT HEIGHT Z = HWIND),
!C* OF POTENTIAL TEMPERATURE DEFICIT (BETWEEN HEIGHT Z = HTEMP AND SURFACE),
!C* OF SPECIFIC HUMIDITY DEFICIT (BETWEEN HEIGHT Z = HTEMP AND SURFACE)
!C* INPUT DATA USED:
!C* ARRAY AR2 (OUTPUT FROM DRAG3 - SUBROUTINE)
!C* AR2(1) - NON-DIMENSIONAL (NORMALIZED BY MONIN-OBUKHOV LENGTH)
!C* HEIGHT OF CONSTANT FLUX LAYER TOP
!C* AR2(5) - ROUGHNESS LENGTH FOR MOMENTUM
!C* AR2(6) - ROUGHNESS LENGTH FOR TEMPERATURE AND HUMIDITY
!C* AR2(8) - INTEGRAL TRANSFER COEFFICIENT FOR MOMENTUM
!C* AR2(9) - INTEGRAL TRANSFER COEFFICIENT FOR TEMPERATURE AND HUMIDITY
!C* ARRAY ARDIN
!C* ARDIN(1) - MODULE OF WIND VELOCITY AT THE TOP OF CONSTANT FLUX LAYER
!C* ARDIN(2) - POTENTIAL TEMPERATURE DEFICIT IN CONSTANT FLUX LAYER
!C* ARDIN(3) - SPECIFIC HUMIDITY DEFICIT IN CONSTANT FLUX LAYER
!C* ARDIN(4) - DIMENSIONAL HEIGHT OF CONSTANT FLUX LAYER TOP
!C* ARDIN(5) - = HWIND
!C* ARDIN(6) - = HTEMP
!C* OUTPUT DATA (ARRAY ARDOUT):
!C* ARDOUT(1) - MODULE OF WIND VELOCITY AT REQUIRED HEIGHT
!C* ARDOUT(2) - POTENTIAL TEMPERATURE at REQUIRED HEIGHT
!C* ARDOUT(3) - SPECIFIC HUMIDITY at REQUIRED HEIGHT
!C* UFWIND - UNIVERSAL FUNCTION FOR WIND VELOCITY
!C* UFTEMP - UNIVERSAL FUNCTION FOR SCALARS
! AR2(1) - NON-DIMENSIONAL CFL HIGHT =
! AR2(2) - RICHARDSON NUMBER =
! AR2(3) - REYNODS NUMBER =
! AR2(4) - LN(ZU/ZT) =
! AR2(5) - DYNAMICAL ROUGHNESS ZU (M) =
! AR2(6) - THERMAL ROUGHNESS ZT (M) =
! AR2(7) - CRITICAL RICHARDSON NUMBER =
! AR2(8) - TRANSFER COEFFICIENT FOR MOMENTUM =
! AR2(9) - TRANSFER COEFFICIENT FR HEAT =
! AR2(10)- COEFFICIENT OF TURBULENCE (KM) AT CFL HIGHT (M**2/S) =
! AR2(11)- ALFT=KT/KM ( KT-COEFFICIENT OF TURBULENCE FOR HEAT) =
subroutine pbldia_new_sheba(AR2,ARDIN,ARDOUT)
use sfx_sheba, only: &
get_psi_sheba => get_psi
real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10)
real,intent(out)::ARDOUT(3) !ARDOUT(4)
! real,parameter::Llim = 50.0 !minimum value of L for stable SL (25)
real,parameter::zetalim = 2. !minimum value of L for stable SL (25)
real,parameter::kappa = 0.4
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.gt.zetalim) zeta=zetalim
! if(L.lt.Llim.and.L.gt.0) L=Llim
L = HIN/zeta ! MO length
call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L)
call get_psi_sheba(psi_m,psi_h,HWIND/L)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
WIND = (USTAR/kappa) * UFWIND
call get_psi_sheba(psi_m,psi_h,HTEMP/L)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
RETURN
END subroutine pbldia_new_sheba
subroutine pbldia_new_most(AR2,ARDIN,ARDOUT)
use sfx_most, only: &
get_psi_most => get_psi
real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10)
real,intent(out)::ARDOUT(3) !ARDOUT(4)
! real,parameter::Llim = 50.0 !minimum value of L for stable SL (25)
real,parameter::zetalim = 2. !minimum value of L for stable SL (25)
real,parameter::kappa = 0.4
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.gt.zetalim) zeta=zetalim
! if(L.lt.Llim.and.L.gt.0) L=Llim
L = HIN/zeta ! MO length
call get_psi_most(psi_m_hs,psi_h_hs,HIN/L)
call get_psi_most(psi_m,psi_h,HWIND/L)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
WIND = (USTAR/kappa) * UFWIND
call get_psi_most(psi_m,psi_h,HTEMP/L)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
! if(L>0)then
! !these functions in very stable stratification give inadequate temperatures at high levels; it's safer to use log profile for T
! UFTEMP = ALOG(HTEMP/HIN)
! endif
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
RETURN
END subroutine pbldia_new_most
subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
!use sfx_esm, only: &
! get_psi_esm => get_psi
real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10)
real,intent(out)::ARDOUT(3) !ARDOUT(4)
! real,parameter::Llim = 50.0 !minimum value of L for stable SL (25)
real,parameter::zetalim = 2. !minimum value of L for stable SL (25)
real,parameter::kappa = 0.4
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.gt.zetalim) zeta=zetalim
! if(L.lt.Llim.and.L.gt.0) L=Llim
L = HIN/zeta ! MO length
call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L)
UFWIND = psi_m
WIND = (USTAR/kappa) * UFWIND
call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
UFTEMP = psi_h
! if(L>0)then
! !these functions in very stable stratification give inadequate temperatures at high levels; it's safer to use log profile for T
! UFTEMP = ALOG(HTEMP/HIN)
! endif
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
RETURN
END subroutine pbldia_new_esm
subroutine get_psi_esm1(psi_m,psi_h,z1,z2,L)
use sfx_esm_param
!< @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions
real, intent(in) :: z1,z2,L !
real an5,D1,gmz1,gtz1,gmz2,gtz2,gmst,gtst,zeta1,zeta2
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)
zeta1 = z1 / L
zeta2 = z2 / L
IF(zeta2 .GE. 0.) THEN
psi_m = alog(z2/z1) + beta_m*(zeta2-zeta1)
psi_h = alog(z2/z1) + beta_m*(zeta2-zeta1)
ELSE
GMZ2 = SQRT(SQRT(1.E0 - alpha_m * zeta2))
GMZ1 = SQRT(SQRT(1.E0 - alpha_m * zeta1))
GMST = SQRT(SQRT(1.E0 - alpha_m * D1))
GTZ2 = SQRT(1.E0 - alpha_h * zeta2)
GTZ1 = SQRT(1.E0 - alpha_h * zeta1)
GTST = SQRT(1.E0 - alpha_h * D1)
IF(zeta2 .GE. D1) THEN
psi_m = FIM(GMZ2) - FIM(GMZ1)
psi_h = FIT(GTZ2) - FIT(GTZ1)
ELSE
psi_m = (3.E0/GMST) * (1.E0 - (D1/zeta2)**(1.E0/3.E0)) + FIM(GMST) - FIM(GMZ1)
psi_h = (3.E0/GTST) * (1.E0 - (D1/zeta2)**(1.E0/3.E0)) + FIT(GTST) - FIT(GTZ1)
ENDIF
ENDIF
end subroutine
REAL FUNCTION FIM(XX)
IMPLICIT NONE
REAL X, XX
X = AMAX1(XX,1.000001E0)
FIM = ALOG((X-1.E0)/(X+1.E0)) + 2.E0*ATAN(X)
RETURN
END function
REAL FUNCTION FIT(XX)
IMPLICIT NONE
REAL X, XX
X = AMAX1(XX,1.000001E0)
FIT = ALOG((X-1.E0)/(X+1.E0))
RETURN
END function
end module pbldia_new_sfx
\ No newline at end of file
diag/sfx_diag_def.fi 0 → 100755
+ 2
− 0
View file @ 393b2460
! sfx diag macro definitions
#define MOST
srcF/sfx_most.f90
+ 1
− 0
View file @ 393b2460
...@@ -23,6 +23,7 @@ module sfx_most ...@@ -23,6 +23,7 @@ module sfx_most
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
public :: get_surface_fluxes public :: get_surface_fluxes
public :: get_surface_fluxes_vec public :: get_surface_fluxes_vec
public :: get_psi
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
......
srcF/sfx_sheba.f90
+ 1
− 0
View file @ 393b2460
...@@ -28,6 +28,7 @@ module sfx_sheba ...@@ -28,6 +28,7 @@ module sfx_sheba
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
public :: get_surface_fluxes public :: get_surface_fluxes
public :: get_surface_fluxes_vec public :: get_surface_fluxes_vec
public :: get_psi
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment