CESE CAM模块学习（2）

今天来尝试解决昨天的问题，最终目标还是找到我的参数化方案的插入位置

cam_out%tbot(i)  = state%t(i,pver)
cam_out%thbot(i) = state%t(i,pver) * state%exner(i,pver)
cam_out%zbot(i)  = state%zm(i,pver)
cam_out%ubot(i)  = state%u(i,pver)
cam_out%vbot(i)  = state%v(i,pver)
cam_out%pbot(i)  = state%pmid(i,pver)

这里的第二个索引，在ppgrid（物理网格化文件）中能找到相关定义

! Grid point resolution parameters

   integer pcols      ! number of columns (max)
   integer psubcols   ! number of sub-columns (max)
   integer pver       ! number of vertical levels
   integer pverp      ! pver + 1

   parameter (pcols     = PCOLS)
   parameter (psubcols  = PSUBCOLS)
   parameter (pver      = PLEV)
   parameter (pverp     = pver + 1  )
!
! start, end indices for chunks owned by a given MPI task
! (set in phys_grid_init).
!
   integer :: begchunk = 0            ! 
   integer :: endchunk = -1           ! 

这里大概就能得出pver真实代表的就是垂直层的最低层索引的结论

接下来继续看module physics_types，查看physics_state相关的定义,因为cam_out来源于physics_state

physics_state下午边开会边看完了 没来得及作记录，但是具有的关键物理变量和一些接口方法都差不多理清了

因为先前的计算方式是参考的风场相关的参数化方案，所以接下来在physpkg文件中查找相关的调用和计算过程

文件中比较关键的两个方法，调用物理过程的两个子例程是

主要以垂直过程为例，查看相关的计算，tphysac子例程的输入参数中有较为熟悉cam_in,cam_out,state,tend,pbuf等变量，囊括了CAM模块所有的重要变量

call t_startf('tphysac_init')
! Associate pointers with physics buffer fields
itim = pbuf_old_tim_idx()


ifld = pbuf_get_index('DTCORE')
call pbuf_get_field(pbuf, ifld, dtcore, start=(/1,1,itim/), kount=(/pcols,pver,1/) )

call pbuf_get_field(pbuf, tini_idx, tini)
call pbuf_get_field(pbuf, qini_idx, qini)
call pbuf_get_field(pbuf, cldliqini_idx, cldliqini)
call pbuf_get_field(pbuf, cldiceini_idx, cldiceini)

ifld = pbuf_get_index('CLD')
call pbuf_get_field(pbuf, ifld, cld, start=(/1,1,itim/),kount=(/pcols,pver,1/))

ifld = pbuf_get_index('AST')
call pbuf_get_field(pbuf, ifld, ast, start=(/1,1,itim/), kount=(/pcols,pver,1/) )

!
! accumulate fluxes into net flux array for spectral dycores
! jrm Include latent heat of fusion for snow
!
do i=1,ncol
   tend%flx_net(i) = tend%flx_net(i) + cam_in%shf(i) + (cam_out%precc(i) &
        + cam_out%precl(i))*latvap*rhoh2o &
        + (cam_out%precsc(i) + cam_out%precsl(i))*latice*rhoh2o
end do

! emission of aerosols at surface
call aerosol_emis_intr (state, cam_in)

if (carma_do_emission) then
   ! carma emissions
   call carma_emission_tend (state, ptend, cam_in, ztodt)
   call physics_update(state, ptend, ztodt, tend)
end if

! get nstep and zero array for energy checker
zero = 0._r8
nstep = get_nstep()
call check_tracers_init(state, tracerint)

! Check if latent heat flux exceeds the total moisture content of the
! lowest model layer, thereby creating negative moisture.

call qneg4('TPHYSAC '       ,lchnk               ,ncol  ,ztodt ,               &
     state%q(1,pver,1),state%rpdel(1,pver) ,cam_in%shf ,         &
     cam_in%lhf , cam_in%cflx )

call t_stopf('tphysac_init')

tphysac子例程的第一个主要部分是进行物理过程的初始化，在调用物理过程之前对一些参数进行准备，主要关于地表通量，气溶胶排放，以及部分函数检查是否有负的湿度

!===================================================
! Source/sink terms for advected tracers.
!===================================================
call t_startf('adv_tracer_src_snk')
! Test tracers

call tracers_timestep_tend(state, ptend, cam_in%cflx, cam_in%landfrac, ztodt)      
call physics_update(state, ptend, ztodt, tend)
call check_tracers_chng(state, tracerint, "tracers_timestep_tend", nstep, ztodt,   &
     cam_in%cflx)

call aoa_tracers_timestep_tend(state, ptend, cam_in%cflx, cam_in%landfrac, ztodt)      
call physics_update(state, ptend, ztodt, tend)
call check_tracers_chng(state, tracerint, "aoa_tracers_timestep_tend", nstep, ztodt,   &
     cam_in%cflx)

! Chemistry calculation
if (chem_is_active()) then
   call chem_timestep_tend(state, ptend, cam_in, cam_out, ztodt, &
        pbuf,  fh2o, fsds)

   call physics_update(state, ptend, ztodt, tend)
   call check_energy_chng(state, tend, "chem", nstep, ztodt, fh2o, zero, zero, zero)
   call check_tracers_chng(state, tracerint, "chem_timestep_tend", nstep, ztodt, &
        cam_in%cflx)
end if
call t_stopf('adv_tracer_src_snk')

第二个部分是关于示踪物质的部分，简单了解不深入

!===================================================
! Vertical diffusion/pbl calculation
! Call vertical diffusion code (pbl, free atmosphere and molecular)
!===================================================

! If CLUBB is called, do not call vertical diffusion, but obukov length and
!   surface friction velocity still need to be computed.  In addition, 
!   surface fluxes need to be updated here for constituents 
if (do_clubb_sgs) then

   call clubb_surface ( state, ptend, ztodt, cam_in, surfric, obklen)
   
   ! Update surface flux constituents 
   call physics_update(state, ptend, ztodt, tend)

else

   call t_startf('vertical_diffusion_tend')
   call vertical_diffusion_tend (ztodt ,state ,cam_in%wsx, cam_in%wsy,   &
        cam_in%shf     ,cam_in%cflx     ,surfric  ,obklen   ,ptend    ,ast    ,&
        cam_in%ocnfrac  , cam_in%landfrac ,        &
        sgh30    ,pbuf )

!------------------------------------------
! Call major diffusion for extended model
!------------------------------------------
if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
   call mspd_intr (ztodt    ,state    ,ptend)
endif

   call physics_update(state, ptend, ztodt, tend)
   call t_stopf ('vertical_diffusion_tend')

endif

这部分是计算垂直扩散过程和边界层高度（PBL）的变化

略过CLUBB模式部分主要逻辑：

• 调用 vertical_diffusion_tend 函数进行垂直扩散的计算，并更新状态变量和时间步趋势。
• 如果是扩展模型（extended model），还会调用主要扩散模块。
• 最后通过 physics_update 函数更新状态变量和时间步趋势

接下来看两个关键一是vertical_diffusion_tend的计算过程二是复习下今天下午看过physics_update