Skip to content
Snippets Groups Projects
Commit fecffcf8 authored by Виктория Суязова's avatar Виктория Суязова Committed by Anna Shestakova
Browse files

added sheba coare, new diag

parent 76f5612f
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 2
− 0
View file @ fecffcf8
......@@ -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
......
diag/pbldia_new_sfx.f90
+ 411
− 310
View file @ fecffcf8
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
srcF/sfx_sheba_noniterative.f90
+ 9
− 30
View file @ fecffcf8
......@@ -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
! --------------------------------------------------------------------------------
......
srcF/sfx_z0m_all_surface.f90
+ 2
− 2
View file @ fecffcf8
......@@ -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
! --------------------------------------------------------------------------------
......
srcF/sfx_z0t_all_surface.f90
+ 2
− 2
View file @ fecffcf8
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment