! 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