diff --git a/CMakeLists.txt b/CMakeLists.txt
index d206ec7567d0e3992f3ffd5677132c4fa9b0c729..39e636a6c21b711c709aa4864de7c52b3f2ccadf 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -93,6 +93,11 @@ set(SOURCES_F
     srcF/sfx_api_inmcm.f90
     srcF/sfx_api_term.f90
 )
+
+set(DIAG 
+    diag/pbldia_new_sfx.f90
+)
+
 set(MAIN_F
         srcF/sfx_main.f90)
 
@@ -162,11 +167,11 @@ endif(USE_CXX OR INCLUDE_CUDA)
 set(SOURCES 
     ${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_F} ${SOURCES_F} ${MAIN_F})
+    ${HEADERS_F} ${SOURCES_F} ${DIAG} ${MAIN_F})
 set(SOURCES_LIB 
     ${MEMPROC_HEADERS_CU} ${MEMPROC_SOURCES_CU} ${MEMPROC_HEADERS_CXX} ${MEMPROC_SOURCES_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 ")
 
diff --git a/diag/compile_simple.sh b/diag/compile_simple.sh
new file mode 100755
index 0000000000000000000000000000000000000000..590a68b45552dfa3949ee04488ef971983ebc2ef
--- /dev/null
+++ b/diag/compile_simple.sh
@@ -0,0 +1,21 @@
+#!/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
diff --git a/diag/diag_sfc_simple.f90 b/diag/diag_sfc_simple.f90
new file mode 100755
index 0000000000000000000000000000000000000000..d59f3833c8013fc233f85696bad4a10d55bd970e
--- /dev/null
+++ b/diag/diag_sfc_simple.f90
@@ -0,0 +1,192 @@
+! 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
diff --git a/diag/pbldia_new_sfx.f90 b/diag/pbldia_new_sfx.f90
new file mode 100755
index 0000000000000000000000000000000000000000..bbcc8cdaa498c4b03e3c4cc29579de2d419831dd
--- /dev/null
+++ b/diag/pbldia_new_sfx.f90
@@ -0,0 +1,252 @@
+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
diff --git a/diag/sfx_diag_def.fi b/diag/sfx_diag_def.fi
new file mode 100755
index 0000000000000000000000000000000000000000..444ddd642c36d9f38a3b50f8ed243e33bb80ce36
--- /dev/null
+++ b/diag/sfx_diag_def.fi
@@ -0,0 +1,2 @@
+! sfx diag macro definitions
+#define MOST
diff --git a/srcF/sfx_most.f90 b/srcF/sfx_most.f90
index 799cf43d671596cd1de134899b53282c23ee2124..f4ce053676793567cd6772a1363c09853d50a45f 100644
--- a/srcF/sfx_most.f90
+++ b/srcF/sfx_most.f90
@@ -23,6 +23,7 @@ module sfx_most
     ! --------------------------------------------------------------------------------
     public :: get_surface_fluxes
     public :: get_surface_fluxes_vec
+    public :: get_psi
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------
diff --git a/srcF/sfx_sheba.f90 b/srcF/sfx_sheba.f90
index c53a396cc768b1e12beca9920a19b4107853ffbb..c33e1a5144fb948ed18a070c835cbdd8319f3f60 100644
--- a/srcF/sfx_sheba.f90
+++ b/srcF/sfx_sheba.f90
@@ -28,6 +28,7 @@ module sfx_sheba
     ! --------------------------------------------------------------------------------
     public :: get_surface_fluxes
     public :: get_surface_fluxes_vec
+    public :: get_psi
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------