每天看一个fortran文件（3） 以及 加入边界层高度到cesm的技术 思路 整理

今天看个关于边界层的

(py27) [chengxl@a3406n08 cam]$ cat pbl_utils.F90 
module 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

一些问题：

要写注释

 ！模块来保存与PBL相关的子程序，这些子程序可以与多个程序一起使用！不同的垂直扩散方案。

!公共的子程序： 计算莫宁奥布霍夫长度

---

 

很多程序都有r8 ,可能是什么定义数据的，我下次也这样写一下。

关于implicit none 和private 就不用说了

---

公共子程序，这些子程序可以被其他的过程调用。

大概有边界层工具初始化程序

计算摩擦速度u_*

计算莫宁奥布霍夫长度L

计算辐射程序

---

这块也没什么好说的，熟悉边界层的人都知道这些含义

重力加速度

冯卡曼常数

干空气定压比热

干空气常数

还有一个和湿空气计算有关的比率

 有意思的是这里，常数后面加个这个玩意

---

 这里是我以前不知道的，interface 

interface 

查了一些资料，可变大小数组用在子程序中需要使用interface

可以参考这个大哥wjfqvi的文章：Fortran的两种显式接口：模块化和interface

写一个subroutine，如果形参给定数组元素个数，可以正常编译运行

subroutine abc(a)
    implicit none 
    integer :: a(3)
    a(1)=1
    a(2)=2 
end subroutine


program main
    implicit none 

    integer :: a(3)=(/5,5,5/)

    call abc(a)

    write(*,*) a 
end 

如果在subroutine里面写数组元素个数，会产生错误

subroutine abc(a)
    implicit none 
    integer :: a(:)
    a(1)=1
    a(2)=2 
end subroutine


program main
    implicit none 

    integer :: a(3)=(/5,5,5/)

    call abc(a)

    write(*,*) a 
end 

错误信息：

$ gfortran 2.f90
2.f90:14:12:

     call abc(a)
            1
Error: Explicit interface required for ‘abc’ at (1): assumed-shape argument


解决方法1：
将过程放在模块中创建显示接口：

module mm

contains 

subroutine abc(a)
    implicit none 
    integer :: a(:)
    a(1)=1
    a(2)=2 
end subroutine
end module

program main
    use mm
    implicit none 

    integer :: a(3)=(/5,5,5/)

    call abc(a)

    write(*,*) a 
end 

解决方法2：
创建接口：

subroutine abc(a)
    implicit none 
    integer :: a(:)
    a(1)=1
    a(2)=2 
end subroutine

program main
    implicit none 

    interface 
        subroutine abc(a)
            implicit none 
            integer :: a(:)
        end subroutine abc 
    end interface

    integer :: a(3)=(/5,5,5/)
    call abc(a)
    write(*,*) a 
end 
————————————————
版权声明：本文为CSDN博主「wjfqvi」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/weixin_43880667/article/details/84836145

解释的好清楚！

---

注意在module 里面也需要先写contains 再写subroutine

---

第一个子程序：就是将一些变量输入进来

---

 

第二个子程序：

和计算摩擦速度的标度有关

我们看到这里也是把xy方向上的摩擦应力给输入进来然后算模

很关键的一个问题是，得看到一个引用这个子程序的程序

这样输入输出很明显，我也要养成这个习惯

---

第三个子程序：计算矢量的摩擦速度

这个程序和上面的差不多，不同的是对于一维数组进行计算

---

 

 

 第四个子程序：计算奥布霍夫长度以及一些通量。边界层的人应该对这些比较熟悉了，反正我是很熟悉。

 计算向量的奥布霍夫长度

 

注意输入文件，很重要 

 

 

---

 

解决问题思路

好，前面也提到这些计算比较清楚，那么一个关键的问题是，输入变量在哪里？这对于我计算边界层高度更为重要。

那么我们看看在哪里call 了这些子程序

我们看见有几个文件都调用了这个函数

  我们看到一个完整的调用是这样的，之后我也就知道这些变量从哪里取了。

 再找下g

 

 很有意思 ，physpkg 是一个和物理过程密切相关的文件

下面我们仿照计算莫宁奥布霍夫长度的计算程序，写一个关于计算边界层高度的程序。

加入代码：

 

 再写一个接口

 

---

 

现在还没有进行调用，先编译一下看看有没有严重的语法错误

最后测试：

 

[chengxl@a3212n16 FHIST_SH_BACK]$ grep -r "pbl_utils" . 
./bld/atm/obj/Depends:clubb_intr.o : clubb_intr.F90 rad_constituents.mod macrop_driver.mod scammod.mod stats_zt_module.mod stats_sfc_module.mod phys_grid.mod cldfrc2m.mod stats_rad_zm_module.mod clubb_precision.mod spmd_utils.mod physics_types.mod stats_zm_module.mod physconst.mod hb_diff.mod ref_pres.mod constituents.mod stats_variables.mod constants_clubb.mod stats_rad_zt_module.modpbl_utils.mod cam_history.mod tropopause.mod co2_cycle.mod zm_conv_intr.mod namelist_utils.mod phys_control.mod time_manager.mod physics_buffer.mod cam_logfile.mod clubb_api_module.mod ppgrid.mod error_code.mod camsrfexch.mod cam_abortutils.mod units.mod   
./bld/atm/obj/Depends:eddy_diff.o : eddy_diff.F90 wv_saturation.mod pbl_utils.mod   
./bld/atm/obj/Depends:eddy_diff_cam.o : eddy_diff_cam.F90 phys_debug_util.mod eddy_diff.mod spmd_utils.mod diffusion_solver.mod physics_types.mod wv_saturation.mod physconst.mod error_messages.mod ref_pres.mod constituents.mod pbl_utils.mod cam_history.mod namelist_utils.mod phys_control.mod physics_buffer.mod time_manager.mod cam_logfile.mod ppgrid.mod coords_1d.mod camsrfexch.mod cam_abortutils.mod units.mod   
./bld/atm/obj/Depends:hb_diff.o : hb_diff.F90 spmd_utils.mod pbl_utils.mod cam_logfile.mod ppgrid.mod   
./bld/atm/obj/Depends:pbl_utils.o : pbl_utils.F90    
./bld/atm/obj/Depends:physpkg.o : physpkg.F90 rad_constituents.mod scammod.mod aoa_tracers.mod convect_shallow.mod phys_grid.mod cam3_ozone_data.mod waccmx_phys_intr.mod iondrag.mod spmd_utils.mod solar_data.mod subcol.mod tracers.mod phys_debug.mod wv_saturation.mod gw_drag.mod convect_deep.mod prescribed_ghg.mod carma_intr.mod conv_water.mod offline_driver.mod nudging.mod flux_avg.mod physconst.mod prescribed_ozone.mod chem_surfvals.mod ref_pres.mod constituents.mod cam_initfiles.mod microp_aero.mod polar_avg.mod cam_grid_support.mod rayleigh_friction.mod pbl_utils.mod prescribed_strataero.mod cam_history.mod subcol_utils.mod radiation.mod string_utils.mod cam_logfile.mod ppgrid.mod chemistry.mod clybry_fam.mod cam_abortutils.mod cam_diagnostics.mod cloud_fraction.mod aero_model.mod macrop_driver.mod check_energy.mod aer_rad_props.mod phys_debug_util.mod ncdio_atm.mod iop_forcing.mod cldfrc2m.mod short_lived_species.mod cam_control_mod.mod mo_lightning.mod physics_types.mod vertical_diffusion.mod dadadj_cam.mod mo_gas_phase_chemdr.mod qneg_module.mod cloud_diagnostics.mod dyn_comp.mod dycore.mod modal_aero_calcsize.mod ghg_data.mod clubb_intr.mod epp_ionization.mod spcam_drivers.mod aircraft_emit.mod mo_apex.mod sslt_rebin.mod aerodep_flx.mod microp_driver.mod prescribed_volcaero.mod phys_gmean.mod charge_neutrality.mod tropopause.mod co2_cycle.mod cam3_aero_data.mod cospsimulator_intr.mod phys_control.mod prescribed_aero.mod time_manager.mod metdata.mod modal_aero_wateruptake.mod physics_buffer.mod radheat.mod qbo.mod rk_stratiform.mod unicon_cam.mod camsrfexch.mod carma_flags_mod.mod   
./bld/atm/obj/Depends:vertical_diffusion.o : vertical_diffusion.F90 rad_constituents.mod spmd_utils.mod wv_saturation.mod physconst.mod hb_diff.mod ref_pres.mod constituents.mod upper_bc.mod pbl_utils.mod cam_history.mod namelist_utils.mod molec_diff.mod cam_logfile.mod ppgrid.mod coords_1d.mod cam_abortutils.mod diffusion_solver.mod physics_types.mod beljaars_drag_cam.mod error_messages.mod trb_mtn_stress_cam.mod eddy_diff_cam.mod co2_cycle.mod phys_control.mod physics_buffer.mod time_manager.mod camsrfexch.mod units.mod   
./bld/atm/obj/Depends:pbl_utils.mod : pbl_utils.o
./bld/atm/obj/Srcfiles:pbl_utils.F90
Binary file ./run/cesm.exe matches

针对这个问题，一位师兄提醒了我有关于程序编译时候，可能依赖关系有没有确定的问题，所以不仅要关注于f90程序，还需要关注到depend文件以及Makefile等等文件里面提示的关系。所以我用grep -r "pbl_util" . 到cime 输出文件里再去找了找依赖关系，很可惜好像不是这个问题，但这个知识对我很有用！感谢师兄

 

 

---

一些问题记录

undefined reference to `calc_pbl_h_vector_'

这种可能是没有use module 造成的,要检查 是否函数拼写错误， 是否public ，是否use ，是否only 。

---

/public/home/chengxl/cesm/components/cam/src/physics/cam/vertical_diffusion.F90(988): error #8284:If the actual argument is scalar, the dummy argument shall be scalar unless the actual argument is of type character or is an element of an array that is not assumed shape, pointer, or polymorphic. 

如果实际参数是标量，则虚拟参数应为标量，除非实际参数是字符类型或不是假定形状、指针或多态的数组元素。

我原本以为重力加速度是变量，没想到是模式里没有计算那么精确，所以我得把我的虚参调成标量。果然在subroutine的申明中我去掉了g_in的数组括号后，编译不再出错。

---

这个问题，是我把subroutine 封装到interface中了，然而实际上我对interface并不熟悉，在经过一晚的调试和讨论后，依然没有解决，最后我回去睡觉了，第二天一早，我就想我为什么非要封装他呢直接用subroutine不就好了？恩，然后就真的解决了这个问题，或者说 我绕开它了。

---

运行中的问题：

我想这应该是由于在某个输出变量表中没有添加我的变量。

有两种解决办法：

第一种移花接木，我放到原先模式本来就要输出的pblh中，这样挂羊头卖狗肉。

首先注释我的输出

 其次把我的PBL_H--->PBLH

使用ncl画图

 

还是存在一点问题，计算的高度偏小，主要原因可能是没有计算fd的变化和N频率的作用。

但是已经成功将方案加入到模式中了！ 

这种移花接木的方法不太好，不推荐！

 

第二种我老老实实把变量表找到，加上这个PBL_H。 关于这个具体的方法我会后面再更新。