今天先来看vertical_diffusion_tend子例程源码，比较经典的参数化方案之一
因为要尝试构建自己的AI参数化方案，对于没有气象学基础的苦命计算机学科学生只能努力先看源码DeBug，再仿写出自己的方案
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 )
在调用时使用了大量cam_in中的变量，包括海洋模块和陆地模块
主要使用变量如下：
ztodt：两个时间步长（2 delta-t）。state：物理状态变量。taux：x方向地表应力（N/m²）。tauy：y方向地表应力（N/m²）。shflx：地表感热通量（W/m²）。cflx：地表组分通量（kg/m²/s）。cldn：新的层积云分数。ocnfrac：海洋分数。landfrac：陆地分数。sgh：地形标准差。ptend：各个参数化趋势。pbuf：物理缓冲区描述符。
其中state基本包括了CAM所有关键物理变量
! ---------------------- ! ! Input-Output Arguments ! ! ---------------------- ! type(physics_ptend), intent(out) :: ptend ! Individual parameterization tendencies type(physics_buffer_desc), pointer :: pbuf(:) ! ---------------- ! ! Output Arguments ! ! ---------------- ! real(r8), intent(out) :: ustar(pcols) ! Surface friction velocity [ m/s ] real(r8), intent(out) :: obklen(pcols) ! Obukhov length [ m ]
输出变量中主要是地表摩擦速度（m/s），Obukhov长度（m），还有ptend梯度
接下来是大量的局部变量定义
logical :: kvinit ! Tell compute_eddy_diff/ caleddy
! to initialize kvh, kvm (uses kvf)
character(128) :: errstring ! Error status for compute_vdiff
real(r8), pointer, dimension(:,:) :: qrl ! LW radiative cooling rate
real(r8), pointer, dimension(:,:) :: wsedl ! Sedimentation velocity
! of stratiform liquid cloud droplet [ m/s ]
integer :: lchnk ! Chunk identifier
integer :: ncol ! Number of atmospheric columns
integer :: i, k, l, m ! column, level, constituent indices
integer :: ierr ! status for allocate/deallocate
real(r8) :: dtk(pcols,pver) ! T tendency from KE dissipation
real(r8), pointer :: tke(:,:) ! Turbulent kinetic energy [ m2/s2 ]
integer(i4),pointer :: turbtype(:,:) ! Turbulent interface types [ no unit ]
real(r8), pointer :: smaw(:,:) ! Normalized Galperin instability function
! ( 0<= <=4.964 and 1 at neutral )
real(r8) :: cgs(pcols,pverp) ! Counter-gradient star [ cg/flux ]
real(r8) :: cgh(pcols,pverp) ! Counter-gradient term for heat
real(r8) :: rztodt ! 1./ztodt [ 1/s ]
real(r8) :: ksrftms(pcols) ! Turbulent mountain stress surface drag coefficient [ kg/s/m2 ]
real(r8) :: tautmsx(pcols) ! U component of turbulent mountain stress [ N/m2 ]
real(r8) :: tautmsy(pcols) ! V component of turbulent mountain stress [ N/m2 ]
real(r8) :: tautotx(pcols) ! U component of total surface stress [ N/m2 ]
real(r8) :: tautoty(pcols) ! V component of total surface stress [ N/m2 ]
real(r8), pointer :: kvh_in(:,:) ! kvh from previous timestep [ m2/s ]
real(r8), pointer :: kvm_in(:,:) ! kvm from previous timestep [ m2/s ]
real(r8), pointer :: kvt(:,:) ! Molecular kinematic conductivity for temperature [ ]
real(r8) :: kvq(pcols,pverp) ! Eddy diffusivity for constituents [ m2/s ]
real(r8) :: kvh(pcols,pverp) ! Eddy diffusivity for heat [ m2/s ]
real(r8) :: kvm(pcols,pverp) ! Eddy diffusivity for momentum [ m2/s ]
real(r8) :: bprod(pcols,pverp) ! Buoyancy production of tke [ m2/s3 ]
real(r8) :: sprod(pcols,pverp) ! Shear production of tke [ m2/s3 ]
real(r8) :: sfi(pcols,pverp) ! Saturation fraction at interfaces [ fraction ]
real(r8) :: sl(pcols,pver)
real(r8) :: qt(pcols,pver)
real(r8) :: slv(pcols,pver)
real(r8) :: sl_prePBL(pcols,pver)
real(r8) :: qt_prePBL(pcols,pver)
real(r8) :: slv_prePBL(pcols,pver)
real(r8) :: slten(pcols,pver)
real(r8) :: qtten(pcols,pver)
real(r8) :: slvten(pcols,pver)
real(r8) :: slflx(pcols,pverp)
real(r8) :: qtflx(pcols,pverp)
real(r8) :: uflx(pcols,pverp)
real(r8) :: vflx(pcols,pverp)
real(r8) :: slflx_cg(pcols,pverp)
real(r8) :: qtflx_cg(pcols,pverp)
real(r8) :: uflx_cg(pcols,pverp)
real(r8) :: vflx_cg(pcols,pverp)
real(r8) :: th(pcols,pver) ! Potential temperature
real(r8) :: topflx(pcols) ! Molecular heat flux at top interface
real(r8) :: wpert(pcols) ! Turbulent wind gusts
real(r8) :: rhoair
real(r8) :: ri(pcols,pver) ! richardson number (HB output)
! for obklen calculation outside HB
real(r8) :: thvs(pcols) ! Virtual potential temperature at surface
real(r8) :: rrho(pcols) ! Reciprocal of density at surface
real(r8) :: khfs(pcols) ! sfc kinematic heat flux [mK/s]
real(r8) :: kqfs(pcols) ! sfc kinematic water vapor flux [m/s]
real(r8) :: kbfs(pcols) ! sfc kinematic buoyancy flux [m^2/s^3]
real(r8) :: ftem(pcols,pver) ! Saturation vapor pressure before PBL
real(r8) :: ftem_prePBL(pcols,pver) ! Saturation vapor pressure before PBL
real(r8) :: ftem_aftPBL(pcols,pver) ! Saturation vapor pressure after PBL
real(r8) :: tem2(pcols,pver) ! Saturation specific humidity and RH
real(r8) :: t_aftPBL(pcols,pver) ! Temperature after PBL diffusion
real(r8) :: tten(pcols,pver) ! Temperature tendency by PBL diffusion
real(r8) :: rhten(pcols,pver) ! RH tendency by PBL diffusion
real(r8) :: qv_aft_PBL(pcols,pver) ! qv after PBL diffusion
real(r8) :: ql_aft_PBL(pcols,pver) ! ql after PBL diffusion
real(r8) :: qi_aft_PBL(pcols,pver) ! qi after PBL diffusion
real(r8) :: s_aft_PBL(pcols,pver) ! s after PBL diffusion
real(r8) :: u_aft_PBL(pcols,pver) ! u after PBL diffusion
real(r8) :: v_aft_PBL(pcols,pver) ! v after PBL diffusion
real(r8) :: qv_pro(pcols,pver)
real(r8) :: ql_pro(pcols,pver)
real(r8) :: qi_pro(pcols,pver)
real(r8) :: s_pro(pcols,pver)
real(r8) :: t_pro(pcols,pver)
real(r8), pointer :: tauresx(:) ! Residual stress to be added in vdiff to correct
real(r8), pointer :: tauresy(:) ! for turb stress mismatch between sfc and atm accumulated.
real(r8) :: ipbl(pcols)
real(r8) :: kpblh(pcols)
real(r8) :: wstarPBL(pcols)
real(r8) :: tpertPBL(pcols)
real(r8) :: qpertPBL(pcols)
real(r8) :: rairi(pcols,pver+1) ! interface gas constant needed for compute_vdiff
real(r8), pointer :: tpert(:)
real(r8), pointer :: qpert(:)
real(r8), pointer :: pblh(:)
real(r8) :: tmp1(pcols) ! Temporary storage
integer :: nstep
real(r8) :: sum1, sum2, sum3, pdelx
real(r8) :: sflx
! Copy state so we can pass to intent(inout) routines that return
! new state instead of a tendency.
real(r8) :: s_tmp(pcols,pver)
real(r8) :: u_tmp(pcols,pver)
real(r8) :: v_tmp(pcols,pver)
real(r8) :: q_tmp(pcols,pver,pcnst)
大量变量都是关于边界层定义的物理变量，等到用到的时候再回来看
! ----------------------- ! ! Main Computation Begins ! ! ----------------------- ! rztodt = 1._r8 / ztodt lchnk = state%lchnk ncol = state%ncol call pbuf_get_field(pbuf, tauresx_idx, tauresx) call pbuf_get_field(pbuf, tauresy_idx, tauresy) call pbuf_get_field(pbuf, tpert_idx, tpert) call pbuf_get_field(pbuf, qpert_idx, qpert) call pbuf_get_field(pbuf, pblh_idx, pblh) call pbuf_get_field(pbuf, turbtype_idx, turbtype)
在计算开始之前，从state和pbuf中获取了一些变量
这里我发现了一个重要的点，大量有用的物理变量其实都存在物理缓存区之中
为了管理pbuf物理缓存，CAM中有大量相关的代码，之后需要留意有哪些需要使用的变量在物理缓存之中
整个vertical_diffusion有大量的参数化方案，因为目的只是学习就只挑eddy_scheme来详细阅读
select case (eddy_scheme)
case ( 'diag_TKE' )
! ---------------------------------------------------------------- !
! At first time step, have eddy_diff.F90:caleddy() use kvh=kvm=kvf !
! This has to be done in compute_eddy_diff after kvf is calculated !
! ---------------------------------------------------------------- !
if( is_first_step() ) then
kvinit = .true.
else
kvinit = .false.
endif
! ---------------------------------------------- !
! Get LW radiative heating out of physics buffer !
! ---------------------------------------------- !
call pbuf_get_field(pbuf, qrl_idx, qrl )
call pbuf_get_field(pbuf, wsedl_idx, wsedl )
! Retrieve eddy diffusivities for heat and momentum from physics buffer
! from last timestep ( if first timestep, has been initialized by inidat.F90 )
call compute_eddy_diff( lchnk , &
pcols , pver , ncol , state%t , state%q(:,:,1) , ztodt , &
state%q(:,:,ixcldliq) , state%q(:,:,ixcldice) , &
state%s , state%rpdel , cldn , qrl , wsedl , &
state%zm , state%zi , state%pmid , state%pint , state%u , state%v , &
taux , tauy , shflx , cflx(:,1) , wstarent , nturb , &
rrho , ustar , pblh , kvm_in , kvh_in , kvm , &
kvh , kvq , cgh , &
cgs , tpert , qpert , wpert , tke , bprod , &
sprod , sfi , kvinit , &
tauresx , tauresy , ksrftms , &
ipbl(:) , kpblh(:) , wstarPBL(:), turbtype , smaw )
代码核心便是调用compute_eddy_diff过程计算热量和动量的涡流扩散率。这些扩散率从上一个时间步获取，并在第一个时间步初始化。
在调用中传递了各种变量，包括模型状态、水汽和液态水的混合比、表面风应力、能量通量等。这个过程会根据给定的参数计算涡流扩散率，并更新这些变量的值。
进入该子例程内部，跳过关于输入输出以及变量准备等具体过程
do iturb = 1, nturb
! Total stress includes 'tms'.
! Here, in computing 'tms', we can use either iteratively changed 'ufd,vfd' or the
! initially given 'u,v' to the PBL scheme. Note that normal stress, 'taux, tauy'
! are not changed by iteration. In order to treat 'tms' in a fully implicit way,
! I am using updated wind, here.
! Compute ustar
call calc_ustar( tfd(:ncol,pver), pmid(:ncol,pver), &
taux(:ncol) - ksrftms(:ncol) * ufd(:ncol,pver), & ! Zonal wind stress
tauy(:ncol) - ksrftms(:ncol) * vfd(:ncol,pver), & ! Meridional wind stress
rrho(:ncol), ustar(:ncol))
minpblh(:ncol) = 100.0_r8 * ustar(:ncol) ! By construction, 'minpblh' is larger than 1 [m] when 'ustar_min = 0.01'.
! Calculate (qt,sl,n2,s2,ri) from a given set of (t,qv,ql,qi,u,v)
call trbintd( &
pcols , pver , ncol , z , ufd , vfd , tfd , pmid , &
s2 , n2 , ri , zi , pi , cldn , qtfd , qvfd , &
qlfd , qi , sfi , sfuh , sflh , slfd , slv , slslope , &
qtslope , chs , chu , cms , cmu )整个方案计算过程最重要部分就是按算涡流扩散率的迭代步骤数nturb循环迭代计算
注释写的挺全，封装也做的挺完整的
首先计算ustar和minpblh
再给定的（tfd，qvfd，u，v）计算（qt，sl，n2，s2，ri）。
! Save initial (i.e., before iterative diffusion) profile of (qt,sl) at each iteration.
! Only necessary for (qt,sl) not (u,v) because (qt,sl) are newly calculated variables.
if( iturb .eq. 1 ) then
qt(:ncol,:) = qtfd(:ncol,:)
sl(:ncol,:) = slfd(:ncol,:)
endif
! Get free atmosphere exchange coefficients. This 'kvf' is not used in UW moist PBL scheme
call austausch_atm( pcols, pver, ncol, ri, s2, kvf )
! Initialize kvh/kvm to send to caleddy, depending on model timestep and iteration number
! This is necessary for 'wstar-based' entrainment closure.
if( iturb .eq. 1 ) then
if( kvinit ) then
! First iteration of first model timestep : Use free tropospheric value or zero.
if( use_kvf ) then
kvh(:ncol,:) = kvf(:ncol,:)
kvm(:ncol,:) = kvf(:ncol,:)
else
kvh(:ncol,:) = 0._r8
kvm(:ncol,:) = 0._r8
endif
else
! First iteration on any model timestep except the first : Use value from previous timestep
kvh(:ncol,:) = kvh_in(:ncol,:)
kvm(:ncol,:) = kvm_in(:ncol,:)
endif
else
! Not the first iteration : Use from previous iteration
kvh(:ncol,:) = kvh_out(:ncol,:)
kvm(:ncol,:) = kvm_out(:ncol,:)
endif
! Calculate eddy diffusivity (kvh_out,kvm_out) and (tke,bprod,sprod) using
! a given (kvh,kvm) which are used only for initializing (bprod,sprod) at
! the first part of caleddy. (bprod,sprod) are fully updated at the end of
! caleddy after calculating (kvh_out,kvm_out)
call caleddy( pcols , pver , ncol , &
slfd , qtfd , qlfd , slv ,ufd , &
vfd , pi , z , zi , &
qflx , shflx , slslope , qtslope , &
chu , chs , cmu , cms ,sfuh , &
sflh , n2 , s2 , ri ,rrho , &
pblh , ustar , &
kvh , kvm , kvh_out , kvm_out , &
tpert , qpert , qrl , kvf , tke , &
wstarent , bprod , sprod , minpblh , wpert , &
tkes , turbtype , sm_aw , &
kbase_o , ktop_o , ncvfin_o , &
kbase_mg , ktop_mg , ncvfin_mg , &
kbase_f , ktop_f , ncvfin_f , &
wet , web , jtbu , jbbu , &
evhc , jt2slv , n2ht , n2hb , &
lwp , opt_depth , radinvfrac, radf , &
wstar , wstar3fact, &
ebrk , wbrk , lbrk , ricl , ghcl , &
shcl , smcl , ghi , shi , smi , &
rii , lengi , wcap , pblhp , cldn , &
ipbl , kpblh , wsedl)
! Calculate errorPBL to check whether PBL produced convergent solutions or not.
if( iturb .eq. nturb ) then
do i = 1, ncol
errorPBL(i) = 0._r8
do k = 1, pver
errorPBL(i) = errorPBL(i) + ( kvh(i,k) - kvh_out(i,k) )**2
end do
errorPBL(i) = sqrt(errorPBL(i)/pver)
end do
end if
! Eddy diffusivities which will be used for the initialization of (bprod,
! sprod) in 'caleddy' at the next iteration step.
if( iturb .gt. 1 .and. iturb .lt. nturb ) then
kvm_out(:ncol,:) = lambda * kvm_out(:ncol,:) + ( 1._r8 - lambda ) * kvm(:ncol,:)
kvh_out(:ncol,:) = lambda * kvh_out(:ncol,:) + ( 1._r8 - lambda ) * kvh(:ncol,:)
endif
! Set nonlocal terms to zero for flux diagnostics, since not used by caleddy.
cgh(:ncol,:) = 0._r8
cgs(:ncol,:) = 0._r8
if( iturb .lt. nturb ) then
! Each time we diffuse the original state
slfd(:ncol,:) = sl(:ncol,:)
qtfd(:ncol,:) = qt(:ncol,:)
ufd(:ncol,:) = u(:ncol,:)
vfd(:ncol,:) = v(:ncol,:)
!------------------------------------------------------------------------
! Check to see if constituent dependent gas constant needed (WACCM-X)
!------------------------------------------------------------------------
if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
rairi(:ncol,1) = rairv(:ncol,1,lchnk)
do k = 2, pver
do i = 1, ncol
rairi(i,k) = 0.5_r8 * (rairv(i,k,lchnk)+rairv(i,k-1,lchnk))
end do
end do
else
rairi(:ncol,:pver+1) = rair
endif
! Diffuse initial profile of each time step using a given (kvh_out,kvm_out)
! In the below 'compute_vdiff', (slfd,qtfd,ufd,vfd) are 'inout' variables.
call compute_vdiff( lchnk , &
pcols , pver , 1 , ncol , pmid , &
pi , rpdel , t , ztodt , taux , &
tauy , shflx , qflx , ntop_turb , nbot_turb , &
kvh_out , kvm_out , kvh_out , cgs , cgh , &
zi , ksrftms , zero , fieldlist_wet, fieldlist_molec, &
ufd , vfd , qtfd , slfd , &
jnk1d , jnk1d , jnk2d , jnk1d , errstring , &
tauresx , tauresy , 0 , cpairv(:,:,lchnk), rairi , .false. )
call handle_errmsg(errstring, subname="compute_vdiff", &
extra_msg="compute_vdiff called from eddy_diff")
! Retrieve (tfd,qvfd,qlfd) from (slfd,qtfd) in order to
! use 'trbintd' at the next iteration.
do k = 1, pver
do i = 1, ncol
! ----------------------------------------------------- !
! Compute the condensate 'qlfd' in the updated profiles !
! ----------------------------------------------------- !
! Option.1 : Assume grid-mean condensate is homogeneously diffused by the moist turbulence scheme.
! This should bs used if 'pseudodiff = .false.' in vertical_diffusion.F90.
! Modification : Need to be check whether below is correct in the presence of ice, qi.
! I should understand why the variation of ice, qi is neglected during diffusion.
templ = ( slfd(i,k) - g*z(i,k) ) / cpair
call qsat( templ, pmid(i,k), es, qs)
ep2 = .622_r8
temps = templ + ( qtfd(i,k) - qs ) / ( cpair / latvap + latvap * qs / ( rair * templ**2 ) )
call qsat( temps, pmid(i,k), es, qs)
qlfd(i,k) = max( qtfd(i,k) - qi(i,k) - qs ,0._r8 )
! Option.2 : Assume condensate is not diffused by the moist turbulence scheme.
! This should bs used if 'pseudodiff = .true.' in vertical_diffusion.F90.
! qlfd(i,k) = ql(i,k)
! ----------------------------- !
! Compute the other 'qvfd, tfd' !
! ----------------------------- !
qvfd(i,k) = max( 0._r8, qtfd(i,k) - qi(i,k) - qlfd(i,k) )
tfd(i,k) = ( slfd(i,k) + latvap * qlfd(i,k) + latsub * qi(i,k) - g*z(i,k)) / cpair
end do
end do
endif
! Debug
! icol = phys_debug_col(lchnk)
! if( icol > 0 .and. get_nstep() .ge. 1 ) then
! write(iulog,*) ' '
! write(iulog,*) 'eddy_diff debug at the end of iteration'
! write(iulog,*) 't, qv, ql, cld, u, v'
! do k = pver-3, pver
! write (iulog,*) k, tfd(icol,k), qvfd(icol,k), qlfd(icol,k), cldn(icol,k), ufd(icol,k), vfd(icol,k)
! end do
! endif
! Debug
end do ! End of 'iturb' iteration如果是第一个迭代，保存初始（即在迭代扩散之前）的（qt，sl）的剖面。
获取自由大气交换系数kvf。
初始化kvh和kvm，以在下一个时间步骤中用于初始化caleddy的（bprod，sprod）。
调用caleddy计算涡流扩散率和（tke，bprod，sprod）。
如果不是最后一次迭代，则更新kvh_out和kvm_out，以在下一个迭代步骤中用于初始化caleddy的（bprod，sprod）。
将非本地项设置为零，以便用于通量诊断。
如果不是最后一次迭代，则将初始状态（slfd，qtfd，ufd，vfd）扩散为下一次迭代的输入。
计算迭代的误差errorPBL
kvq(:ncol,:) = kvh_out(:ncol,:)
! Compute 'wstar' within the PBL for use in the future convection scheme.
do i = 1, ncol
if( ipbl(i) .eq. 1._r8 ) then
wstarPBL(i) = max( 0._r8, wstar(i,1) )
else
wstarPBL(i) = 0._r8
endif
end do
! --------------------------------------------------------------- !
! Writing for detailed diagnostic analysis of UW moist PBL scheme !
! --------------------------------------------------------------- !
call outfld( 'UW_errorPBL', errorPBL, pcols, lchnk )
call outfld( 'UW_n2', n2, pcols, lchnk )
call outfld( 'UW_s2', s2, pcols, lchnk )
call outfld( 'UW_ri', ri, pcols, lchnk )
call outfld( 'UW_sfuh', sfuh, pcols, lchnk )
call outfld( 'UW_sflh', sflh, pcols, lchnk )
call outfld( 'UW_sfi', sfi, pcols, lchnk )
call outfld( 'UW_cldn', cldn, pcols, lchnk )
call outfld( 'UW_qrl', qrl, pcols, lchnk )
call outfld( 'UW_ql', qlfd, pcols, lchnk )
call outfld( 'UW_chu', chu, pcols, lchnk )
call outfld( 'UW_chs', chs, pcols, lchnk )
call outfld( 'UW_cmu', cmu, pcols, lchnk )
call outfld( 'UW_cms', cms, pcols, lchnk )
call outfld( 'UW_tke', tke, pcols, lchnk )
call outfld( 'UW_wcap', wcap, pcols, lchnk )
call outfld( 'UW_bprod', bprod, pcols, lchnk )
call outfld( 'UW_sprod', sprod, pcols, lchnk )
call outfld( 'UW_kvh', kvh_out, pcols, lchnk )
call outfld( 'UW_kvm', kvm_out, pcols, lchnk )
call outfld( 'UW_pblh', pblh, pcols, lchnk )
call outfld( 'UW_pblhp', pblhp, pcols, lchnk )
call outfld( 'UW_tpert', tpert, pcols, lchnk )
call outfld( 'UW_qpert', qpert, pcols, lchnk )
call outfld( 'UW_wpert', wpert, pcols, lchnk )
call outfld( 'UW_ustar', ustar, pcols, lchnk )
call outfld( 'UW_tkes', tkes, pcols, lchnk )
call outfld( 'UW_minpblh', minpblh, pcols, lchnk )
call outfld( 'UW_turbtype', real(turbtype,r8), pcols, lchnk )
call outfld( 'UW_kbase_o', kbase_o, pcols, lchnk )
call outfld( 'UW_ktop_o', ktop_o, pcols, lchnk )
call outfld( 'UW_ncvfin_o', ncvfin_o, pcols, lchnk )
call outfld( 'UW_kbase_mg', kbase_mg, pcols, lchnk )
call outfld( 'UW_ktop_mg', ktop_mg, pcols, lchnk )
call outfld( 'UW_ncvfin_mg', ncvfin_mg, pcols, lchnk )
call outfld( 'UW_kbase_f', kbase_f, pcols, lchnk )
call outfld( 'UW_ktop_f', ktop_f, pcols, lchnk )
call outfld( 'UW_ncvfin_f', ncvfin_f, pcols, lchnk )
call outfld( 'UW_wet', wet, pcols, lchnk )
call outfld( 'UW_web', web, pcols, lchnk )
call outfld( 'UW_jtbu', jtbu, pcols, lchnk )
call outfld( 'UW_jbbu', jbbu, pcols, lchnk )
call outfld( 'UW_evhc', evhc, pcols, lchnk )
call outfld( 'UW_jt2slv', jt2slv, pcols, lchnk )
call outfld( 'UW_n2ht', n2ht, pcols, lchnk )
call outfld( 'UW_n2hb', n2hb, pcols, lchnk )
call outfld( 'UW_lwp', lwp, pcols, lchnk )
call outfld( 'UW_optdepth', opt_depth, pcols, lchnk )
call outfld( 'UW_radfrac', radinvfrac, pcols, lchnk )
call outfld( 'UW_radf', radf, pcols, lchnk )
call outfld( 'UW_wstar', wstar, pcols, lchnk )
call outfld( 'UW_wstar3fact', wstar3fact, pcols, lchnk )
call outfld( 'UW_ebrk', ebrk, pcols, lchnk )
call outfld( 'UW_wbrk', wbrk, pcols, lchnk )
call outfld( 'UW_lbrk', lbrk, pcols, lchnk )
call outfld( 'UW_ricl', ricl, pcols, lchnk )
call outfld( 'UW_ghcl', ghcl, pcols, lchnk )
call outfld( 'UW_shcl', shcl, pcols, lchnk )
call outfld( 'UW_smcl', smcl, pcols, lchnk )
call outfld( 'UW_gh', ghi, pcols, lchnk )
call outfld( 'UW_sh', shi, pcols, lchnk )
call outfld( 'UW_sm', smi, pcols, lchnk )
call outfld( 'UW_ria', rii, pcols, lchnk )
call outfld( 'UW_leng', lengi, pcols, lchnk )
return
end subroutine compute_eddy_diff最后是该方案的结尾工作
- 将一个数组的值赋给另一个数组：
kvq(:ncol,:) = kvh_out(:ncol,:) - 计算并保存在 PBL 内部的
wstar值以备将来对对流方案的使用 - 输出一系列诊断字段，包括但不限于
errorPBL、n2、s2、ri、sfuh、sflh、sfi、cldn、qrl、qlfd、chu、chs、cmu、cms、tke、wcap、bprod、sprod、kvh_out、kvm_out、pblh、pblhp、tpert、qpert、wpert、ustar、tkes、minpblh、turbtype、kbase_o、kbase_mg、kbase_f等。这些输出用于进一步分析 UW 湿润边界层方案的性能和行为。
该方案与我的研究内容无关，但也是第一次较为完整的阅读完了整个过程，在不深究物理过程的细节下，整个代码对一个外行来说阅读难度不高，所以希望自己打起精神来，要做的东西很难，但是总有突破口。凡事都是有代价的，越难的东西做出来就越有价值。
接下来想看看深对流过程的代码，再和学长学姐讨论下
扫一扫
2089
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
