diff --git a/diag/compile_simple.sh b/diag/compile_simple.sh
index 590a68b45552dfa3949ee04488ef971983ebc2ef..c4db5e0fa3d10aa046e1dbe977614de1dc38f660 100755
--- a/diag/compile_simple.sh
+++ b/diag/compile_simple.sh
@@ -2,6 +2,7 @@
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_sheba_noit_param.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_most_param.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_log_param.f90
@@ -10,6 +11,7 @@ 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_sheba_noniterative.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_most.f90
gfortran -c -cpp -Wuninitialized ../srcF/sfx_log.f90
@@ -17,5 +19,5 @@ 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
+gfortran -o a.exe sfx_phys_const.o sfx_sheba_param.o sfx_sheba_noit_param.o sfx_esm_param.o sfx_most_param.o sfx_log_param.o sfx_data.o sfx_surface.o sfx_sheba.o sfx_sheba_noniterative.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
index d59f3833c8013fc233f85696bad4a10d55bd970e..89acfd47b1e7e94d79e3f7eda9d5f0699494387b 100755
--- a/diag/diag_sfc_simple.f90
+++ b/diag/diag_sfc_simple.f90
@@ -8,6 +8,9 @@ use sfx_data
use sfx_sheba, only: &
get_surface_fluxes_vec_sheba => get_surface_fluxes_vec, &
numericsType_sheba => numericsType
+use sfx_sheba_noniterative, only: &
+ get_surface_fluxes_vec_sheba_noit => get_surface_fluxes_vec, &
+ numericsType_sheba_noit => numericsType
use sfx_esm, only: &
get_surface_fluxes_vec_esm => get_surface_fluxes_vec, &
numericsType_esm => numericsType
@@ -26,6 +29,9 @@ use sfx_phys_const
type(meteoDataVecType) :: meteo !< meteorological data (input)
type(meteoDataType) :: meteo_cell
type(sfxDataVecType) :: sfx !< surface fluxes (output)
+#ifdef SHEBA_NOIT
+type(numericsType_sheba_noit) :: numerics_sheba_noit !< surface flux module (SHEBA) numerics parameters
+#endif
#ifdef SHEBA
type(numericsType_sheba) :: numerics_sheba !< surface flux module (SHEBA) numerics parameters
#endif
@@ -54,14 +60,14 @@ integer j,k
! user-defined input parameters for wind, temperature and humidity interpolation in the surface layer
!mandatory:
-u = 5.0
+u = 10.
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
+h = 50. ! input height
zw = 10 ! height to which wind interpolation is performed
zt = 2 ! height to wich temperature and humidity interpolation is performed
@@ -133,6 +139,9 @@ meteo%z0_m = z0m
#ifdef SHEBA
call get_surface_fluxes_vec_sheba(sfx, meteo, numerics_sheba, 1)
#endif
+#ifdef SHEBA_NOIT
+call get_surface_fluxes_vec_sheba_noit(sfx, meteo, numerics_sheba_noit, 1)
+#endif
#ifdef ESM
call get_surface_fluxes_vec_esm(sfx, meteo, numerics_esm, 1)
#endif
@@ -152,13 +161,14 @@ ar2(9) = sfx%Ct(1)
else
ar2(1) = h / L
-ar2(5) = z0m
-ar2(6) = z0h
+!ar2(5) = z0m
+!ar2(6) = z0h
ar2(8) = cm
ar2(9) = ct
endif
+print *, sfx%zeta,-1*1005.*sfx%Cm(1)*sfx%Ct(1)*u*dt*1.3
ardin(1) = u
ardin(2) = dt
@@ -175,6 +185,9 @@ ardin(6) = zt
#ifdef SHEBA
call pbldia_new_sheba(ar2(:),ardin(:),ardout(:))
#endif
+#ifdef SHEBA_NOIT
+call pbldia_new_sheba_noit(ar2(:),ardin(:),ardout(:))
+#endif
#ifdef MOST
call pbldia_new_most(ar2(:),ardin(:),ardout(:))
#endif
diff --git a/diag/pbldia_new_sfx.f90 b/diag/pbldia_new_sfx.f90
index 684d89dc8ac78008865a0a0d0f2bcda5523968a9..e11e3b51b66605bea31c930d52d35349c49b30c7 100755
--- a/diag/pbldia_new_sfx.f90
+++ b/diag/pbldia_new_sfx.f90
@@ -84,13 +84,78 @@ use sfx_sheba, only: &
return
end subroutine pbldia_new_sheba
+subroutine pbldia_new_sheba_noit(AR2,ARDIN,ARDOUT)
+use sfx_sheba_noniterative, only: &
+ get_psi_stable_sheba => get_psi_stable
+
+ real,intent(in):: AR2(11),ARDIN(6)
+ real,intent(out)::ARDOUT(3)
+ real,parameter::zetalim = 2. !maximum value of z/L for stable SL
+ 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.ne.0)then
+
+ if(zeta.gt.zetalim) zeta=zetalim
+ L = HIN/zeta
+
+ if(zeta.gt.0)then ! stable stratification
+
+
+ call get_psi_stable_sheba(psi_m_hs,psi_h_hs,HIN/L,HIN/L)
+ call get_psi_stable_sheba(psi_m,psi_h,HWIND/L,HWIND/L)
+ UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
+ call get_psi_stable_sheba(psi_m,psi_h,HTEMP/L,HTEMP/L)
+ UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
+
+ else ! unstable stratification (functions from ESM)
+
+ call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L)
+ UFWIND = psi_m
+ call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
+ UFTEMP = psi_h
+
+ endif
+
+ else
+
+! neutral stratification (z/L=0)
+
+ UFWIND = ALOG(HWIND/HIN)
+ UFTEMP = ALOG(HTEMP/HIN)
+
+ endif
+
+ WIND = (USTAR/kappa) * UFWIND
+ 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_noit
+
subroutine pbldia_new_most(AR2,ARDIN,ARDOUT)
use sfx_most, only: &
get_psi_most => get_psi
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
- real,parameter::zetalim = 1. !maximum value of z/L for stable SL
+ real,parameter::zetalim = 2. !maximum value of z/L for stable SL
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
@@ -137,7 +202,7 @@ use sfx_most, only: &
subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
- real,parameter::zetalim = 1. !maximum value of z/L for stable SL
+ real,parameter::zetalim = 2. !maximum value of z/L for stable SL
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
diff --git a/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90
index 41e0dae3dcef8f568f803b969d5770685ba11d71..178d89ef707274e185ebb6f601233e69581d0e93 100644
--- a/srcF/sfx_sheba_noniterative.f90
+++ b/srcF/sfx_sheba_noniterative.f90
@@ -27,6 +27,7 @@ module sfx_sheba_noniterative
! --------------------------------------------------------------------------------
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
+ public :: get_psi_stable
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
@@ -38,16 +39,21 @@ module sfx_sheba_noniterative
#if defined(INCLUDE_CXX)
type, BIND(C), public :: sfx_sheba_noit_param_C
- real(C_FLOAT) :: kappa
+ real(C_FLOAT) :: kappa
real(C_FLOAT) :: Pr_t_0_inv
real(C_FLOAT) :: Pr_t_inf_inv
- real(C_FLOAT) :: alpha_m
+ real(C_FLOAT) :: alpha_m
real(C_FLOAT) :: alpha_h
real(C_FLOAT) :: alpha_h_fix
real(C_FLOAT) :: Rib_max
real(C_FLOAT) :: gamma
real(C_FLOAT) :: zeta_a
+ real(C_FLOAT) :: a_m
+ real(C_FLOAT) :: b_m
+ real(C_FLOAT) :: a_h
+ real(C_FLOAT) :: b_h
+ real(C_FLOAT) :: c_h
end type
type, BIND(C), public :: sfx_sheba_noit_numericsType_C
@@ -70,21 +76,6 @@ module sfx_sheba_noniterative
type(sfx_phys_constants) :: constants
END SUBROUTINE c_sheba_noit_compute_flux
- SUBROUTINE c_sheba_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, &
- name="c_sheba_compute_flux")
- use sfx_data
- use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR
- Import :: sfx_sheba_param_C, sfx_sheba_numericsType_C
- implicit none
- INTEGER(C_INT) :: grid_size
- type(C_PTR), value :: sfx
- type(C_PTR), value :: meteo
- type(sfx_sheba_param_C) :: model_param
- type(sfx_surface_param) :: surface_param
- type(sfx_sheba_numericsType_C) :: numerics
- type(sfx_phys_constants) :: constants
- END SUBROUTINE c_sheba_compute_flux
-
END INTERFACE
#endif
@@ -104,23 +95,12 @@ contains
sfx_model_param%Rib_max = Rib_max
sfx_model_param%gamma = gamma
sfx_model_param%zeta_a = zeta_a
- end subroutine set_c_struct_sfx_sheba_noit_param_values
-#endif
-
-#if defined(INCLUDE_CXX)
- subroutine set_c_struct_sfx_sheba_param_values(sfx_model_param)
- type (sfx_sheba_param_C), intent(inout) :: sfx_model_param
- sfx_model_param%kappa = kappa
- sfx_model_param%Pr_t_0_inv = Pr_t_0_inv
-
- sfx_model_param%alpha_m = alpha_m
- sfx_model_param%alpha_h = alpha_h
sfx_model_param%a_m = a_m
sfx_model_param%b_m = b_m
sfx_model_param%a_h = a_h
sfx_model_param%b_h = b_h
sfx_model_param%c_h = c_h
- end subroutine set_c_struct_sfx_sheba_param_values
+ end subroutine set_c_struct_sfx_sheba_noit_param_values
#endif
! --------------------------------------------------------------------------------
@@ -150,12 +130,6 @@ contains
type (sfx_phys_constants) :: phys_constants
numerics_c%maxiters_convection = numerics%maxiters_convection
-
- type (sfx_sheba_param_C) :: model_param
- type (sfx_surface_param) :: surface_param
- type (sfx_sheba_numericsType_C) :: numerics_c
- type (sfx_phys_constants) :: phys_constants
-
numerics_c%maxiters_charnock = numerics%maxiters_charnock
phys_constants%Pr_m = Pr_m;
@@ -168,7 +142,7 @@ contains
call set_sfx_vec_c(sfx, sfx_c)
meteo_c_ptr = C_LOC(meteo_c)
sfx_c_ptr = C_LOC(sfx_c)
-
+
call c_sheba_noit_compute_flux(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
#else
do i = 1, n
@@ -245,8 +219,6 @@ contains
real fval !< just a shortcut for partial calculations
- real :: C1,A1,A2,lne,lnet,Ribl
-
#ifdef SFX_CHECK_NAN
real NaN
@@ -316,21 +288,9 @@ contains
! --- restrict bulk Ri value
! *: note that this value is written to output
- Rib = min(Rib, Rib_max)
-
- Ribl = (Rib*Pr_t_0_inv) * (1 - z0_t / h) / ((1 - z0_m / h)**2)
-
- call get_psi_stable(psi_m, psi_h, zeta_a, zeta_a)
- call get_psi_stable(psi0_m, psi0_h, zeta_a * z0_m / h, zeta_a * z0_t / h)
-
- lne = log(h/z0_m)
- lnet = log(h/z0_t)
- C1 = (lne**2)/lnet
- A1 = ((lne - psi_m + psi0_m)**(2*(gamma-1))) &
-& / ((zeta_a**(gamma-1))*((lnet-(psi_h-psi0_h)*Pr_t_0_inv)**(gamma-1)))
- A2 = ((lne - psi_m + psi0_m)**2) / (lnet-(psi_h-psi0_h)*Pr_t_0_inv) - C1
+! Rib = min(Rib, Rib_max)
- zeta = C1 * Ribl + A1 * A2 * (Ribl**gamma)
+ call get_zeta(zeta, Rib, h, z0_m, z0_t)
call get_psi_stable(psi_m, psi_h, zeta, zeta)
call get_psi_stable(psi0_m, psi0_h, zeta * z0_m / h, zeta * z0_t / h)
@@ -397,6 +357,29 @@ contains
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
+ subroutine get_zeta(zeta, Rib, h, z0_m, z0_t)
+ real,intent(out) :: zeta
+ real,intent(in) :: Rib, h, z0_m, z0_t
+
+ real :: Ribl, C1, A1, A2, lne, lnet
+ real :: psi_m, psi_h, psi0_m, psi0_h
+
+ Ribl = (Rib * Pr_t_0_inv) * (1 - z0_t / h) / ((1 - z0_m / h)**2)
+
+ call get_psi_stable(psi_m, psi_h, zeta_a, zeta_a)
+ call get_psi_stable(psi0_m, psi0_h, zeta_a * z0_m / h, zeta_a * z0_t / h)
+
+ lne = log(h/z0_m)
+ lnet = log(h/z0_t)
+ C1 = (lne**2)/lnet
+ A1 = ((lne - psi_m + psi0_m)**(2*(gamma-1))) &
+& / ((zeta_a**(gamma-1))*((lnet-(psi_h-psi0_h)*Pr_t_0_inv)**(gamma-1)))
+ A2 = ((lne - psi_m + psi0_m)**2) / (lnet-(psi_h-psi0_h)*Pr_t_0_inv) - C1
+
+ zeta = C1 * Ribl + A1 * A2 * (Ribl**gamma)
+
+ end subroutine get_zeta
+
! convection universal functions shortcuts
! --------------------------------------------------------------------------------
function f_m_conv(zeta)
@@ -469,7 +452,7 @@ contains
log((2.0 * x_h + c_h - q_h) / (2.0 * x_h + c_h + q_h)) - &
log((c_h - q_h) / (c_h + q_h)))
- end subroutine
+ end subroutine get_psi_stable
subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)