问题：把大象放到冰箱需要几步
三步。
打开冰箱；
放入大象；
关上冰箱；
把公式放到模式里要几步：
三步。
打开模式；
放入公式；
关上模式；
第一步 打开模式相关代码文件
就是找到计算边界层高度的程序：
这怎么知道呀
@login04 cam]$ cd /public/home/chengxl/cesm/components/cam
@login04 cam]$ find -name '*pbl*'
./src/physics/cam/pbl_utils.F90找到一个文件
[chengxl@login04 cam]$ vim ./src/physics/cam/pbl_utils.F90module pbl_utils
!-----------------------------------------------------------------------!
! Module to hold PBL-related subprograms that may be used with multiple !
! different vertical diffusion schemes. !
! !
! Public subroutines: !
!
! calc_obklen !
! !
!------------------ History --------------------------------------------!
! Created: Apr. 2012, by S. Santos !
!-----------------------------------------------------------------------!
use shr_kind_mod, only: r8 => shr_kind_r8
implicit none
private
! Public Procedures
!----------------------------------------------------------------------!
! Excepting the initialization procedure, these are elemental
! procedures, so they can accept scalars or any dimension of array as
! arguments, as long as all arguments have the same number of
! elements.
public pbl_utils_init
public calc_ustar
public calc_obklen
public virtem
public compute_radf
public austausch_atm
real(r8), parameter :: ustar_min = 0.01_r8
real(r8) :: g ! acceleration of gravity
real(r8) :: vk ! Von Karman's constant
real(r8) :: cpair ! specific heat of dry air
real(r8) :: rair ! gas constant for dry air
real(r8) :: zvir ! rh2o/rair - 1
!------------------------------------------------------------------------!
! Purpose: Compilers aren't creating optimized vector versions of !
! elemental routines, so we'll explicitly create them and bind !
! them via an interface for transparent use !
!------------------------------------------------------------------------!
interface calc_ustar
module procedure calc_ustar_scalar
module procedure calc_ustar_vector
end interface
interface calc_obklen
module procedure calc_obklen_scalar
module procedure calc_obklen_vector
end interface
interface virtem
module procedure virtem_vector1D
module procedure virtem_vector2D ! Used in hb_diff.F90
end interface
contains
subroutine pbl_utils_init(g_in,vk_in,cpair_in,rair_in,zvir_in)
!-----------------------------------------------------------------------!
! Purpose: Set constants to be used in calls to later functions !
!-----------------------------------------------------------------------!
real(r8), intent(in) :: g_in ! acceleration of gravity
real(r8), intent(in) :: vk_in ! Von Karman's constant
real(r8), intent(in) :: cpair_in ! specific heat of dry air
real(r8), intent(in) :: rair_in ! gas constant for dry air
real(r8), intent(in) :: zvir_in ! rh2o/rair - 1
g = g_in
vk = vk_in
cpair = cpair_in
rair = rair_in
zvir = zvir_in
end subroutine pbl_utils_init
subroutine calc_ustar_scalar( t, pmid, taux, tauy, &
rrho, ustar)
!-----------------------------------------------------------------------!
! Purpose: Calculate ustar and bottom level density (necessary for !
! Obukhov length calculation). !
!-----------------------------------------------------------------------!
real(r8), intent(in) :: t ! surface temperature
real(r8), intent(in) :: pmid ! midpoint pressure (bottom level)
real(r8), intent(in) :: taux ! surface u stress [N/m2]
real(r8), intent(in) :: tauy ! surface v stress [N/m2]
real(r8), intent(out) :: rrho ! 1./bottom level density
real(r8), intent(out) :: ustar ! surface friction velocity [m/s]
rrho = rair * t / pmid
ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min )
end subroutine calc_ustar_scalar
subroutine calc_ustar_vector(n, t, pmid, taux, tauy, &
rrho, ustar)
!-----------------------------------------------------------------------!
! Purpose: Calculate ustar and bottom level density (necessary for !
! Obukhov length calculation). !
!-----------------------------------------------------------------------!
integer, intent(in) :: n ! Length of vectors
real(r8), intent(in) :: t(n) ! surface temperature
real(r8), intent(in) :: pmid(n) ! midpoint pressure (bottom level)
real(r8), intent(in) :: taux(n) ! surface u stress [N/m2]
real(r8), intent(in) :: tauy(n) ! surface v stress [N/m2]
real(r8), intent(out) :: rrho(n) ! 1./bottom level density
real(r8), intent(out) :: ustar(n) ! surface friction velocity [m/s]
rrho = rair * t / pmid
ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min )
end subroutine calc_ustar_vector
subroutine calc_obklen_scalar( ths, thvs, qflx, shflx, rrho, ustar, &
khfs, kqfs, kbfs, obklen)
!-----------------------------------------------------------------------!
! Purpose: Calculate Obukhov length and kinematic fluxes. !
!-----------------------------------------------------------------------!
real(r8), intent(in) :: ths ! potential temperature at surface [K]
real(r8), intent(in) :: thvs ! virtual potential temperature at surface
real(r8), intent(in) :: qflx ! water vapor flux (kg/m2/s)
real(r8), intent(in) :: shflx ! surface heat flux (W/m2)
real(r8), intent(in) :: rrho ! 1./bottom level density [ m3/kg ]
real(r8), intent(in) :: ustar ! Surface friction velocity [ m/s ]
real(r8), intent(out) :: khfs ! sfc kinematic heat flux [mK/s]
real(r8), intent(out) :: kqfs ! sfc kinematic water vapor flux [m/s]
real(r8), intent(out) :: kbfs ! sfc kinematic buoyancy flux [m^2/s^3]
real(r8), intent(out) :: obklen ! Obukhov length
! Need kinematic fluxes for Obukhov:
khfs = shflx*rrho/cpair
kqfs = qflx*rrho
kbfs = khfs + zvir*ths*kqfs
! Compute Obukhov length:
obklen = -thvs * ustar**3 / (g*vk*(kbfs + sign(1.e-10_r8,kbfs)))
end subroutine calc_obklen_scalar
subroutine calc_obklen_vector(n, ths, thvs, qflx, shflx, rrho, ustar, &
khfs, kqfs, kbfs, obklen)
!-----------------------------------------------------------------------!
! Purpose: Calculate Obukhov length and kinematic fluxes. !
!-----------------------------------------------------------------------!
integer, intent(in) :: n ! Length of vectors
real(r8), intent(in) :: ths(n) ! potential temperature at surface [K]
real(r8), intent(in) :: thvs(n) ! virtual potential temperature at surface
real(r8), intent(in) :: qflx(n) ! water vapor flux (kg/m2/s)
real(r8), intent(in) :: shflx(n) ! surface heat flux (W/m2)
real(r8), intent(in) :: rrho(n) ! 1./bottom level density [ m3/kg ]
real(r8), intent(in) :: ustar(n) ! Surface friction velocity [ m/s ]
real(r8), intent(out) :: khfs(n) ! sfc kinematic heat flux [mK/s]
real(r8), intent(out) :: kqfs(n) ! sfc kinematic water vapor flux [m/s]
real(r8), intent(out) :: kbfs(n) ! sfc kinematic buoyancy flux [m^2/s^3]
real(r8), intent(out) :: obklen(n) ! Obukhov length
! Need kinematic fluxes for Obukhov:
khfs = shflx*rrho/cpair
kqfs = qflx*rrho
kbfs = khfs + zvir*ths*kqfs
! Compute Obukhov length:
obklen = -thvs * ustar**3 / (g*vk*(kbfs + sign(1.e-10_r8,kbfs)))
end subroutine calc_obklen_vector
subroutine virtem_vector1D(n, t,q, virtem)
!-----------------------------------------------------------------------!
! Purpose: Calculate virtual temperature from temperature and specific !
! humidity. !
!-----------------------------------------------------------------------!
integer, intent(in) :: n ! vector length
real(r8), intent(in) :: t(n), q(n)
real(r8), intent(out):: virtem(n)
virtem = t * (1.0_r8 + zvir*q)
end subroutine virtem_vector1D
subroutine virtem_vector2D(n, m, t, q, virtem)
!-----------------------------------------------------------------------!
! Purpose: Calculate virtual temperature from temperature and specific !
! humidity. !
!-----------------------------------------------------------------------!
integer, intent(in) :: n, m ! vector lengths
real(r8), intent(in) :: t(n,m), q(n,m)
real(r8), intent(out):: virtem(n,m)
virtem = t * (1.0_r8 + zvir*q)
end subroutine virtem_vector2D
subroutine compute_radf( choice_radf, i, pcols, pver, ncvmax, ncvfin, ktop, qmin, &
ql, pi, qrlw, g, cldeff, zi, chs, lwp_CL, opt_depth_CL, &
radinvfrac_CL, radf_CL )
! -------------------------------------------------------------------------- !
! Purpose: !
! Calculate cloud-top radiative cooling contribution to buoyancy production. !
! Here, 'radf' [m2/s3] is additional buoyancy flux at the CL top interface !
! associated with cloud-top LW cooling being mainly concentrated near the CL !
! top interface ( just below CL top interface ). Contribution of SW heating !
! within the cloud is not included in this radiative buoyancy production !
! since SW heating is more broadly distributed throughout the CL top layer. !
! -------------------------------------------------------------------------- !
!-----------------!
! Input variables !
!-----------------!
character(len=6), intent(in) :: choice_radf ! Method for calculating radf
integer, intent(in) :: i ! Index of current column
integer, intent(in) :: pcols ! Number of atmospheric columns
integer, intent(in) :: pver ! Number of atmospheric layers
integer, intent(in) :: ncvmax ! Max numbers of CLs (perhaps equal to pver)
integer, intent(in) :: ncvfin(pcols) ! Total number of CL in column
integer, intent(in) :: ktop(pcols, ncvmax) ! ktop for current column
real(r8), intent(in) :: qmin ! Minimum grid-mean LWC counted as clouds [kg/kg]
real(r8), intent(in) :: ql(pcols, pver) ! Liquid water specific humidity [ kg/kg ]
real(r8), intent(in) :: pi(pcols, pver+1) ! Interface pressures [ Pa ]
real(r8), intent(in) :: qrlw(pcols, pver) ! Input grid-mean LW heating rate : [ K/s ] * cpair * dp = [ W/kg*Pa ]
real(r8), intent(in) :: g ! Gravitational acceleration
real(r8), intent(in) :: cldeff(pcols,pver) ! Effective Cloud Fraction [fraction]
real(r8), intent(in) :: zi(pcols, pver+1) ! Interface heights [ m ]
real(r8), intent(in) :: chs(pcols, pver+1) ! Buoyancy coeffi. saturated sl (heat) coef. at all interfaces.
!------------------!
! Output variables !
!------------------!
real(r8), intent(out) :: lwp_CL(ncvmax) ! LWP in the CL top layer [ kg/m2 ]
real(r8), intent(out) :: opt_depth_CL(ncvmax) ! Optical depth of the CL top layer
real(r8), intent(out) :: radinvfrac_CL(ncvmax) ! Fraction of LW radiative cooling confined in the top portion of CL
real(r8), intent(out) :: radf_CL(ncvmax) ! Buoyancy production at the CL top due to radiative cooling [ m2/s3 ]
!-----------------!
! Local variables !
!-----------------!
integer :: kt, ncv
real(r8) :: lwp, opt_depth, radinvfrac, radf
!-----------------!
! Begin main code !
!-----------------!
lwp_CL = 0._r8
opt_depth_CL = 0._r8
radinvfrac_CL = 0._r8
radf_CL = 0._r8
! ---------------------------------------- !
! Perform do loop for individual CL regime !
! ---------------------------------------- !
do ncv = 1, ncvfin(i)
kt = ktop(i,ncv)
!-----------------------------------------------------!
! Compute radf for each CL regime and for each column !
!-----------------------------------------------------!
if( choice_radf .eq. 'orig' ) then
if( ql(i,kt) .gt. qmin .and. ql(i,kt-1) .lt. qmin ) then
lwp = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
! Approximate LW cooling fraction concentrated at the inversion by using
! polynomial approx to exact formula 1-2/opt_depth+2/(exp(opt_depth)-1))
radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
radf = qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) ! Cp*radiative cooling = [ W/kg ]
radf = max( radinvfrac * radf * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 ) * chs(i,kt)
! We can disable cloud LW cooling contribution to turbulence by uncommenting:
! radf = 0._r8
end if
elseif( choice_radf .eq. 'ramp' ) then
lwp = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
radinvfrac = max(cldeff(i,kt)-cldeff(i,kt-1),0._r8) * radinvfrac
radf = qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) ! Cp*radiative cooling [W/kg]
radf = max( radinvfrac * radf * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 ) * chs(i,kt)
elseif( choice_radf .eq. 'maxi' ) then
! Radiative flux divergence both in 'kt' and 'kt-1' layers are included
! 1. From 'kt' layer
lwp = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
radf = max( radinvfrac * qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 )
! 2. From 'kt-1' layer and add the contribution from 'kt' layer
lwp = ql(i,kt-1) * ( pi(i,kt) - pi(i,kt-1) ) / g
opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth) + opt_depth**2 )
radf = radf + max( radinvfrac * qrlw(i,kt-1) / ( pi(i,kt-1) - pi(i,kt) ) * ( zi(i,kt-1) - zi(i,kt) ), 0._r8 )
radf = max( radf, 0._r8 ) * chs(i,kt)
endif
lwp_CL(ncv) = lwp
opt_depth_CL(ncv) = opt_depth
radinvfrac_CL(ncv) = radinvfrac
radf_CL(ncv) = radf
end do ! ncv = 1, ncvfin(i)
end subroutine compute_radf
subroutine austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf)
!---------------------------------------------------------------------- !
! !
! Purpose: Computes exchange coefficients for free turbulent flows. !
! !
! Method: !
! !
! The free atmosphere diffusivities are based on standard mixing length !
! forms for the neutral diffusivity multiplied by functns of Richardson !
! number. K = l^2 * |dV/dz| * f(Ri). The same functions are used for !
! momentum, potential temperature, and constitutents. !
! !
! The stable Richardson num function (Ri>0) is taken from Holtslag and !
! Beljaars (1989), ECMWF proceedings. f = 1 / (1 + 10*Ri*(1 + 8*Ri)) !
! The unstable Richardson number function (Ri<0) is taken from CCM1. !
! f = sqrt(1 - 18*Ri) !
! !
! Author: B. Stevens (rewrite, August 2000) !
! !
!---------------------------------------------------------------------- !
! --------------- !
! Input arguments !
! --------------- !
integer, intent(in) :: pcols ! Atmospheric columns dimension size
integer, intent(in) :: ncol ! Number of atmospheric columns
integer, intent(in) :: pver ! Number of atmospheric layers
integer, intent(in) :: ntop ! Top layer for calculation
integer, intent(in) :: nbot ! Bottom layer for calculation
real(r8), intent(in) :: ml2(pver+1) ! Mixing lengths squared
real(r8), intent(in) :: s2(pcols,pver) ! Shear squared
real(r8), intent(in) :: ri(pcols,pver) ! Richardson no
! ---------------- !
! Output arguments !
! ---------------- !
real(r8), intent(out) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers
! --------------- !
! Local Variables !
! --------------- !
real(r8) :: fofri ! f(ri)
real(r8) :: kvn ! Neutral Kv
integer :: i ! Longitude index
integer :: k ! Vertical index
real(r8), parameter :: zkmin = 0.01_r8 ! Minimum kneutral*f(ri).
! ----------------------- !
! Main Computation Begins !
! ----------------------- !
kvf(:ncol,:) = 0.0_r8
! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm.
do k = ntop, nbot - 1
do i = 1, ncol
if( ri(i,k) < 0.0_r8 ) then
fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) )
else
fofri = 1.0_r8 / ( 1.0_r8 + 10.0_r8 * ri(i,k) * ( 1.0_r8 + 8.0_r8 * ri(i,k) ) )
end if
kvn = ml2(k) * sqrt(s2(i,k))
kvf(i,k+1) = max( zkmin, kvn * fofri )
end do
end do
end subroutine austausch_atm
end module pbl_utils第二步 加入代码
主要方法就是这个公式了
program jisuanbianjiecenggaodu
implicit none
real :: t_c=27
real :: tau=0.01
real :: hsb=10
print*, " height of PBL "
print*,"==================================="
call height_of_PBL(t_c,tau,hsb)
end program jisuanbianjiecenggaodu
subroutine height_of_PBL(T_c,tau_star,hsb_star)
!--------------------------------------!
!purpose:calculate the height of PBL !
!mothod from the zili(2005) paper !
! !
!--------------------------------------!
real :: C_R =0.6
real :: C_CN =1.36
real :: C_NS =0.51
real :: C_N = 0.1
real :: C_f = 1.0
real :: k = 0.4
real :: pi = 3.1415926
real :: lat =10
real :: N
real :: f
real :: varphi
real :: omega
real :: g = 9.8
real :: T_0
real :: T_c
real :: beta_b
real :: tau_star
real :: hsb_star
real :: h_1
real :: h_2
real :: h_3
real :: h_123
real :: h
omega = 2*pi/360*4.167*0.001
varphi = lat*2*pi/360
f = 2*omega*sin(varphi)
N = 0 !brunt-vaisala pinglv
g =9.8
T_0 = 273+T_c
beta_b = g/T_0
h_1 = (f**2)/((C_R**2)*tau_star)
h_2 = (N*abs(f))/(tau_star*(C_CN**2))
h_3 = abs(f*hsb_star*beta_b)/(C_CN**2*tau_star**2)
h_123 = h_1+h_2+h_3
h = sqrt(1/(h_1+h_2+h_3))
print*,"pi = ",pi,"k=",k,"C_f=",C_f,"omega=",omega,"h = ",h
end subroutine 第三步 运行看看
@server02 ~]$ gfortran height.f90
@server02 ~]$ ./a.out
height of PBL
=========================================================================
pi = 3.14159250 ,k = 0.400000006 ,C_f = 1.00000000
omega = 7.27278675E-05 ,h = 4.73462009 
1164
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
