diff --git a/CMakeLists.txt b/CMakeLists.txt
index d15265549cf8abac34f5d1b5b1ac5f9f06402b18..98b57e977e9cdd3321525f2586d2f1923702d3f8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -95,6 +95,8 @@ set(SOURCES_F
srcF/sfx_sheba_param.f90
srcF/sfx_sheba_noniterative.f90
srcF/sfx_sheba_noit_param.f90
+ srcF/sfx_sheba_coare.f90
+ srcF/sfx_sheba_coare_param.f90
srcF/sfx_fc_wrapper.F90
srcF/sfx_api_inmcm.f90
srcF/sfx_api_term.f90
diff --git a/diag/pbldia_new_sfx.f90 b/diag/pbldia_new_sfx.f90
index d30f498f3e08366ae0e8d495a3b1d4e43255fb29..390eaad182fea5c76bb72bf21d5b0b3653524baf 100755
--- a/diag/pbldia_new_sfx.f90
+++ b/diag/pbldia_new_sfx.f90
@@ -1,313 +1,414 @@
module pbldia_new_sfx
-implicit none
-private
-
-public :: pbldia_new_sheba,pbldia_new_sheba_noit,&
-pbldia_new_most,pbldia_new_esm
-
-real,parameter,private::kappa=0.4,R=287.,appa=0.287
-
-contains
-
-! in this module 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(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) - deficit of POTENTIAL TEMPERATURE between REQUIRED HEIGHT and the surface
-!C* ARDOUT(3) - deficit of SPECIFIC HUMIDITY between REQUIRED HEIGHT and the surface
-!C* UFWIND - UNIVERSAL FUNCTION FOR WIND VELOCITY
-!C* UFTEMP - UNIVERSAL FUNCTION FOR SCALARS
-
-
-subroutine pbldia_new_sheba(AR2,ARDIN,ARDOUT)
-use sfx_sheba, only: &
- get_psi_sheba => get_psi
-
- 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
-
- 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)
- call get_psi_sheba(psi_m,psi_h,HTEMP/L)
- UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
-
- else
-
- 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
-
-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 = 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
-
- 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)
- call get_psi_most(psi_m,psi_h,HTEMP/L)
- UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
-
- else
-
- 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_most
-
-subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
- 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
-
- 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
-
- else
-
- 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_esm
-
-
-subroutine get_psi_esm1(psi_m,psi_h,z1,z2,L)
-
-use sfx_esm_param
-
-real, intent(out) :: psi_m, psi_h
-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 get_psi_esm1
-
-! functions fim,fit were copied from the original pbldia
-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
+ implicit none
+ private
+
+ public :: pbldia_new_sheba,pbldia_new_sheba_noit,&
+ pbldia_new_most,pbldia_new_esm,pbldia_new_log,pbldia_new_sheba_coare
+
+ real,parameter,private::kappa=0.4,R=287.,appa=0.287
+
+ contains
+
+ ! in this module 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(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) - deficit of POTENTIAL TEMPERATURE between REQUIRED HEIGHT and the surface
+ !C* ARDOUT(3) - deficit of SPECIFIC HUMIDITY between REQUIRED HEIGHT and the surface
+ !C* UFWIND - UNIVERSAL FUNCTION FOR WIND VELOCITY
+ !C* UFTEMP - UNIVERSAL FUNCTION FOR SCALARS
+
+
+ subroutine pbldia_new_log(AR2,ARDIN,ARDOUT,lat)
+ real,intent(in):: AR2(11),ARDIN(6),lat
+ real,intent(out)::ARDOUT(3)
+ real hwind,htemp,ustar,tstar,qstar
+ real dteta,dq,wind,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)
+
+ WIND = (USTAR/kappa)*ALOG(hwind/HIN)
+ DTETA = (TSTAR/kappa)*ALOG(htemp/HIN)
+ DQ = (QSTAR/kappa)*ALOG(htemp/HIN)
+
+
+ ARDOUT(1) = ARDIN(1)+WIND
+ ARDOUT(2) = DTETA+ardin(2)
+ ARDOUT(3) = DQ+ardin(3)
+
+
+ return
+ end subroutine pbldia_new_log
+
+ subroutine pbldia_new_sheba(AR2,ARDIN,ARDOUT,lat)
+ use sfx_sheba, only: get_psi_sheba => get_psi
+
+ real,intent(in):: AR2(11),ARDIN(6),lat
+ real,intent(out)::ARDOUT(3)
+ real,parameter::zetalim = 1. !maximum value of z/L for
+ ! stable SL for wind
+ 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
+ real Rib,ksi,wg,f,L_lim
+
+ HWIND = ARDIN(5)
+ HTEMP = ARDIN(6)
+ HIN = ardin(4)
+ zeta = AR2(1)
+ Rib = AR2(2)
+
+ USTAR = AR2(8) * ARDIN(1)
+ TSTAR = AR2(9) * ARDIN(2)
+ QSTAR = AR2(9) * ARDIN(3)
+
+ if(zeta.ne.0)then
+
+
+ if(zeta.gt.zetalim) L_lim = HIN/zetalim
+
+ L = HIN/zeta
+
+ call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L_lim)
+ call get_psi_sheba(psi_m,psi_h,HWIND/L_lim)
+ UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
+ call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L)
+ call get_psi_sheba(psi_m,psi_h,HTEMP/L)
+ UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
+
+ else
+
+ 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
+
+ subroutine pbldia_new_sheba_coare(AR2,ARDIN,ARDOUT,lat)
+ use sfx_sheba_coare, only: &
+ get_psi_coare => get_psi_a, &
+ get_psi_stable_sheba => get_psi_stable
+
+ real,intent(in):: AR2(11),ARDIN(6),lat
+ real,intent(out)::ARDOUT(3)
+ real,parameter::zetalim = 1. !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,L_lim
+
+ 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) L_lim = HIN/zetalim
+ L = HIN/zeta
+
+ if(zeta.gt.0)then ! stable stratification
+
+
+ call get_psi_stable_sheba(psi_m_hs,psi_h_hs,HIN/L_lim,HIN/L)
+ call get_psi_stable_sheba(psi_m,psi_h,HWIND/L_lim,HTEMP/L)
+
+ else
+
+ call get_psi_coare(psi_m_hs,psi_h_hs,HIN/L,HIN/L)
+ call get_psi_coare(psi_m,psi_h,HWIND/L,HTEMP/L)
+
+ endif
+
+ UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
+ UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
+
+ else
+
+ 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_coare
+
+
+ subroutine pbldia_new_sheba_noit(AR2,ARDIN,ARDOUT,lat)
+ use sfx_sheba_noniterative, only: &
+ get_psi_stable_sheba => get_psi_stable
+
+ real,intent(in):: AR2(11),ARDIN(6),lat
+ real,intent(out)::ARDOUT(3)
+ real,parameter::zetalim = 1. !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
+ real Rib,wg,ksi,f, L_lim
+
+ HWIND = ARDIN(5)
+ HTEMP = ARDIN(6)
+ HIN = ardin(4)
+ zeta = AR2(1)
+ Rib = AR2(2)
+
+ USTAR = AR2(8) * ARDIN(1)
+ TSTAR = AR2(9) * ARDIN(2)
+ QSTAR = AR2(9) * ARDIN(3)
+
+ if(zeta.ne.0)then
+
+ if(zeta.gt.zetalim) L_lim=HIN/zetalim
+ L = HIN/zeta
+
+ if(zeta.gt.0)then ! stable stratification
+
+
+ call get_psi_stable_sheba(psi_m_hs,psi_h_hs,HIN/L_lim,HIN/L)
+ call get_psi_stable_sheba(psi_m,psi_h,HWIND/L_lim,HTEMP/L)
+ UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
+ 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,lat)
+ use sfx_most, only: &
+ get_psi_most => get_psi
+
+ real,intent(in):: AR2(11),ARDIN(6),lat
+ real,intent(out)::ARDOUT(3)
+ real,parameter::zetalim = 1. !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
+ real Rib,ksi,wg,f,L_lim
+
+ HWIND = ARDIN(5)
+ HTEMP = ARDIN(6)
+ HIN = ardin(4)
+ zeta = AR2(1)
+ Rib = AR2(2)
+
+ USTAR = AR2(8) * ARDIN(1)
+ TSTAR = AR2(9) * ARDIN(2)
+ QSTAR = AR2(9) * ARDIN(3)
+
+ if(zeta.ne.0)then
+
+ if(zeta.gt.zetalim) L_lim=HIN/zetalim
+ L = HIN/zeta
+
+ call get_psi_most(psi_m_hs,psi_h_hs,HIN/L_lim)
+ call get_psi_most(psi_m,psi_h,HWIND/L_lim)
+ UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
+ call get_psi_most(psi_m_hs,psi_h_hs,HIN/L)
+ call get_psi_most(psi_m,psi_h,HTEMP/L)
+ UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
+
+ else
+
+ 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_most
+
+ subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT,lat)
+ real,intent(in):: AR2(11),ARDIN(6),lat
+ real,intent(out)::ARDOUT(3)
+ real,parameter::zetalim = 1. !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_lim,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) L_lim=HIN/zetalim
+ L = HIN/zeta
+
+ call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L_lim)
+ UFWIND = psi_m
+ call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
+ UFTEMP = psi_h
+
+ else
+
+ 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_esm
+
+
+ subroutine get_psi_esm1(psi_m,psi_h,z1,z2,L)
+
+ use sfx_esm_param
+
+ real, intent(out) :: psi_m, psi_h
+ 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 get_psi_esm1
+
+ ! functions fim,fit were copied from the original pbldia
+ 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/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90
index 36fce06ca2c9acae0dc256a6b58b27759f040fe8..503d1ea61599f9e6ca601c8df4c158116f607674 100644
--- a/srcF/sfx_sheba_noniterative.f90
+++ b/srcF/sfx_sheba_noniterative.f90
@@ -288,32 +288,11 @@ contains
call get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, &
h0_m, h0_t, B)
- ! --- define snow consentration
- !call get_sigma(sigma_r, sigma_w, rho_air, rho_s)
- !call get_w_snow(w_snow, sigma_w, g, d_s, nu_air)
- !call get_h_salt(h_salt, u_dyn0)
- !call get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)
- !call get_S_mean(S_mean, S_salt, h_salt, h, w_snow, u_dyn0)
-
- !deltaS=S_salt-S_mean
-
- !Ri_sn = (g * sigma_r * deltaS * h) / U**2
- !Ri_sn=0.0
! --- 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)
- !do i = 1, numerics%maxiters_convection
- !if (surface_type == surface_snow) then
- !write(*,*) 'RIsnow', Rib, Ri_sn
- ! Rib=Rib+Ri_sn
- !endif
-
+
call get_zeta(zeta, Rib, h, z0_m, z0_t)
!write(*,*) 'get_zeta', zeta, h
@@ -407,15 +386,15 @@ contains
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)
- !write(*,*) 'Smean_0ut', S_mean
sfx = sfxDataType(zeta = zeta, Rib = Rib, &
- Re = Linv, B = B, z0_m = z0_m, z0_t = z0_t, &
- Rib_conv_lim = S_mean, &
- Cm = Cm, Ct = Ct, Km = S_salt, Pr_t_inv = Udyn)
+ 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)
+ !write(*,*) 'Smean_0ut', S_mean
+ !sfx = sfxDataType(zeta = zeta, Rib = Rib, &
+ ! Re = Linv, B = B, z0_m = z0_m, z0_t = z0_t, &
+ ! Rib_conv_lim = S_mean, &
+ ! Cm = Cm, Ct = Ct, Km = S_salt, Pr_t_inv = Udyn)
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
diff --git a/srcF/sfx_z0m_all_surface.f90 b/srcF/sfx_z0m_all_surface.f90
index 249ce699569617e050b827cfec42e57894756620..8f2f11ec734ceec8b89fc30f261605cc99825448 100644
--- a/srcF/sfx_z0m_all_surface.f90
+++ b/srcF/sfx_z0m_all_surface.f90
@@ -84,7 +84,7 @@ module sfx_z0m_all_surface
! --- define dynamic velocity in neutral conditions
u_dyn0 = Uc / c
- write(*,*) 'ch', z0_m, u_dyn0
+ !!write(*,*) 'ch', z0_m, u_dyn0
end subroutine
! --------------------------------------------------------------------------------
@@ -204,7 +204,7 @@ subroutine get_dynamic_roughness_map(z0_m, u_dyn0, U, h, z0m_map)
z0_m=z0m_map
h0_m = h / z0_m
u_dyn0 = U * kappa / log(h0_m)
- write(*,*) 'map', z0_m, u_dyn0
+ !write(*,*) 'map', z0_m, u_dyn0
end subroutine
! --------------------------------------------------------------------------------
diff --git a/srcF/sfx_z0t_all_surface.f90 b/srcF/sfx_z0t_all_surface.f90
index c6581f3b82f83bed460de75bfed918af5a7be352..503eeeb5e47b2159d2445541136022bac32e89e2 100644
--- a/srcF/sfx_z0t_all_surface.f90
+++ b/srcF/sfx_z0t_all_surface.f90
@@ -49,7 +49,7 @@ module sfx_z0t_all_surface
B = min(B, B_max_land)
z0_t = z0_m / exp(B)
- write(*,*) 'kl_land', z0_t, B
+ !write(*,*) 'kl_land', z0_t, B
end subroutine
! --------------------------------------------------------------------------------
@@ -75,7 +75,7 @@ module sfx_z0t_all_surface
B = min(B, B_max_ocean)
z0_t = z0_m / exp(B)
- write(*,*) 'kl_water', z0_t, B
+ !(*,*) 'kl_water', z0_t, B
end subroutine