module eddy_diff_cam use shr_kind_mod, only: i4 => shr_kind_i4, r8 => shr_kind_r8 use ppgrid, only: pcols, pver, pverp use cam_logfile, only: iulog use cam_abortutils, only: endrun use physconst, only: gravit, cpair, rair, zvir, latvap, latice, karman use diffusion_solver, only: vdiff_selector use eddy_diff, only: ncvmax use time_manager, only: is_first_step use physics_buffer, only: physics_buffer_desc use spmd_utils, only: masterproc use phys_control, only: phys_getopts implicit none private public :: eddy_diff_readnl public :: eddy_diff_register public :: eddy_diff_init public :: eddy_diff_tend ! Is UNICON switched on (and thus interacting with eddy_diff via pbuf)? logical :: unicon_is_on ! Number of iterations for solution integer, parameter :: nturb = 5 ! Logical switches for moist mixing ratio diffusion type(vdiff_selector) :: fieldlist_wet ! Logical switches for molecular diffusion ! (Molecular diffusion is not done here.) type(vdiff_selector) :: fieldlist_molec integer :: ntop_eddy, nbot_eddy ! Cloud mass constituent indices integer :: ixcldliq, ixcldice ! input pbuf field indices integer :: qrl_idx = -1 integer :: wsedl_idx = -1 ! output pbuf field indices for UNICON integer :: bprod_idx = -1 integer :: ipbl_idx = -1 integer :: kpblh_idx = -1 integer :: wstarPBL_idx = -1 integer :: tkes_idx = -1 integer :: went_idx = -1 ! Mixing lengths squared. ! Used for computing free air diffusivity. real(r8) :: ml2(pver+1) ! Various namelist options to limit or tweak the effects of eddy diffusion. ! Pressure defining the bottom of the upper atmosphere for kvh scaling (Pa) real(r8) :: kv_top_pressure = 0._r8 ! Eddy diffusivity scale factor for upper atmosphere real(r8) :: kv_top_scale = 1._r8 ! Eddy diffusivity scale factor for the free troposphere real(r8) :: kv_freetrop_scale = 1._r8 ! The following all have to be set in all cases. real(r8), parameter :: unset_r8 = huge(1._r8) ! Maximum master length for diag_TKE real(r8) :: eddy_lbulk_max = unset_r8 ! Maximum dissipation length for diag_TKE real(r8) :: eddy_leng_max = unset_r8 ! Bottom pressure level (hPa) for eddy_leng_max real(r8) :: eddy_max_bot_pressure = unset_r8 ! Moist entrainment enhancement param real(r8) :: eddy_moist_entrain_a2l = unset_r8 contains subroutine eddy_diff_readnl(nlfile) use namelist_utils, only: find_group_name use units, only: getunit, freeunit use spmd_utils, only: masterprocid, mpi_real8, mpicom use shr_log_mod, only: errMsg => shr_log_errMsg ! filepath for file containing namelist input character(len=*), intent(in) :: nlfile ! file unit and error code integer :: unitn, ierr character(len=*), parameter :: subname = 'eddy_diff_readnl' namelist /eddy_diff_nl/ kv_top_pressure, kv_top_scale, & kv_freetrop_scale, eddy_lbulk_max, eddy_leng_max, & eddy_max_bot_pressure, eddy_moist_entrain_a2l if (masterproc) then unitn = getunit() open( unitn, file=trim(nlfile), status='old' ) call find_group_name(unitn, 'eddy_diff_nl', status=ierr) if (ierr == 0) then read(unitn, eddy_diff_nl, iostat=ierr) end if if (ierr /= 0) then call endrun(subname // ':: ERROR reading namelist') end if close(unitn) call freeunit(unitn) end if call mpi_bcast(kv_top_pressure, 1, mpi_real8, masterprocid, mpicom, ierr) if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error") call mpi_bcast(kv_top_scale, 1, mpi_real8, masterprocid, mpicom, ierr) if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error") call mpi_bcast(kv_freetrop_scale, 1, mpi_real8, masterprocid, mpicom, ierr) if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error") call mpi_bcast(eddy_lbulk_max, 1, mpi_real8, masterprocid, mpicom, ierr) if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error") call mpi_bcast(eddy_leng_max, 1, mpi_real8, masterprocid, mpicom, ierr) if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error") call mpi_bcast(eddy_max_bot_pressure, 1, mpi_real8, masterprocid, mpicom, ierr) if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error") call mpi_bcast(eddy_moist_entrain_a2l, 1, mpi_real8, masterprocid, mpicom, ierr) if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error") end subroutine eddy_diff_readnl subroutine eddy_diff_register() use physics_buffer, only: pbuf_add_field, dtype_r8, dtype_i4 character(len=16) :: shallow_scheme ! Check for UNICON and add relevant pbuf entries. call phys_getopts(shallow_scheme_out=shallow_scheme) unicon_is_on = (shallow_scheme == "UNICON") if (unicon_is_on) then call pbuf_add_field('bprod', 'global', dtype_r8, (/pcols,pverp/), bprod_idx) call pbuf_add_field('ipbl', 'global', dtype_i4, (/pcols/), ipbl_idx) call pbuf_add_field('kpblh', 'global', dtype_i4, (/pcols/), kpblh_idx) call pbuf_add_field('wstarPBL', 'global', dtype_r8, (/pcols/), wstarPBL_idx) call pbuf_add_field('tkes', 'global', dtype_r8, (/pcols/), tkes_idx) call pbuf_add_field('went', 'global', dtype_r8, (/pcols/), went_idx) end if end subroutine eddy_diff_register subroutine eddy_diff_init(pbuf2d, ntop_eddy_in, nbot_eddy_in) use error_messages, only: handle_errmsg use cam_history, only: addfld, add_default, horiz_only use constituents, only: cnst_get_ind use ref_pres, only: pref_mid use diffusion_solver, only: new_fieldlist_vdiff, vdiff_select use eddy_diff, only: init_eddy_diff use physics_buffer, only: pbuf_set_field, pbuf_get_index type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! Physics buffer integer, intent(in) :: ntop_eddy_in ! Top interface level to which eddy vertical diffusivity is applied ( = 1 ) integer, intent(in) :: nbot_eddy_in ! Bottom interface level to which eddy vertical diffusivity is applied ( = pver ) character(len=128) :: errstring real(r8) :: leng_max(pver) integer :: k logical :: history_amwg ntop_eddy = ntop_eddy_in nbot_eddy = nbot_eddy_in do k = 1, pver if (pref_mid(k) <= eddy_max_bot_pressure*1.e2_r8) then leng_max(k) = eddy_leng_max else leng_max(k) = 40.e3_r8 end if end do if (masterproc) then write(iulog,*)'init_eddy_diff: nturb=',nturb write(iulog,*)'init_eddy_diff: eddy_leng_max=',eddy_leng_max,' lbulk_max=',eddy_lbulk_max do k = 1,pver write(iulog,*)'init_eddy_diff:',k,pref_mid(k),'leng_max=',leng_max(k) end do end if call init_eddy_diff(pver, gravit, cpair, rair, zvir, & latvap, latice, ntop_eddy, nbot_eddy, karman, & eddy_lbulk_max, leng_max, & eddy_moist_entrain_a2l, errstring) call handle_errmsg(errstring, subname="init_eddy_diff") ! Set the square of the mixing lengths. ml2(1:ntop_eddy) = 0._r8 do k = ntop_eddy + 1, nbot_eddy ml2(k) = 30.0_r8**2 end do ml2(nbot_eddy+1:pver+1) = 0._r8 ! Get fieldlists to pass to diffusion solver. fieldlist_wet = new_fieldlist_vdiff(1) fieldlist_molec = new_fieldlist_vdiff(1) call handle_errmsg(vdiff_select(fieldlist_wet,'s'), & subname="vdiff_select") call handle_errmsg(vdiff_select(fieldlist_wet,'q',1), & subname="vdiff_select") call handle_errmsg(vdiff_select(fieldlist_wet,'u'), & subname="vdiff_select") call handle_errmsg(vdiff_select(fieldlist_wet,'v'), & subname="vdiff_select") ! Cloud mass constituents call cnst_get_ind('CLDLIQ', ixcldliq) call cnst_get_ind('CLDICE', ixcldice) ! Input pbuf fields qrl_idx = pbuf_get_index('QRL') wsedl_idx = pbuf_get_index('WSEDL') ! Initialize output pbuf fields if (is_first_step() .and. unicon_is_on) then call pbuf_set_field(pbuf2d, bprod_idx, 1.0e-5_r8) call pbuf_set_field(pbuf2d, ipbl_idx, 0 ) call pbuf_set_field(pbuf2d, kpblh_idx, 1 ) call pbuf_set_field(pbuf2d, wstarPBL_idx, 0.0_r8) call pbuf_set_field(pbuf2d, tkes_idx, 0.0_r8) call pbuf_set_field(pbuf2d, went_idx, 0.0_r8) end if ! Scheme-specific default output. call phys_getopts(history_amwg_out=history_amwg) call addfld('WGUSTD', horiz_only, 'A', 'm/s', 'wind gusts from turbulence' ) if (history_amwg) then call add_default( 'WGUSTD ', 1, ' ' ) end if ! ------------------------------------------------------------------- ! ! Writing outputs for detailed analysis of UW moist turbulence scheme ! ! ------------------------------------------------------------------- ! call addfld( 'BPROD', ['ilev'], 'A', 'm2/s3', 'Buoyancy Production' ) call addfld( 'SFI', ['ilev'], 'A', '1', 'Interface-layer sat frac' ) call addfld( 'SPROD', ['ilev'], 'A', 'm2/s3', 'Shear Production' ) call addfld('UW_errorPBL',horiz_only,'A', 'm2/s', 'Error function of UW PBL') call addfld('UW_n2', ['lev'], 'A', 's-2', 'Buoyancy Frequency, LI') call addfld('UW_s2', ['lev'], 'A', 's-2', 'Shear Frequency, LI') call addfld('UW_ri', ['lev'], 'A', '1', 'Interface Richardson Number, I') call addfld('UW_sfuh', ['lev'], 'A', '1', 'Upper-Half Saturation Fraction, L') call addfld('UW_sflh', ['lev'], 'A', '1', 'Lower-Half Saturation Fraction, L') call addfld('UW_sfi', ['ilev'], 'A', '1', 'Interface Saturation Fraction, I') call addfld('UW_cldn', ['lev'], 'A', '1', 'Cloud Fraction, L') call addfld('UW_qrl', ['lev'], 'A', 'gravity W/m2', 'LW cooling rate, L') call addfld('UW_ql', ['lev'], 'A', 'kg/kg', 'ql(LWC), L') call addfld('UW_chu', ['ilev'], 'A', 'gravity kg/J', 'Buoyancy Coefficient, chu, I') call addfld('UW_chs', ['ilev'], 'A', 'gravity kg/J', 'Buoyancy Coefficient, chs, I') call addfld('UW_cmu', ['ilev'], 'A','gravity/kg/kg', 'Buoyancy Coefficient, cmu, I') call addfld('UW_cms', ['ilev'], 'A','gravity/kg/kg', 'Buoyancy Coefficient, cms, I') call addfld('UW_tke', ['ilev'], 'A', 'm2/s2', 'TKE, I') call addfld('UW_wcap', ['ilev'], 'A', 'm2/s2', 'Wcap, I') call addfld('UW_bprod', ['ilev'], 'A', 'm2/s3', 'Buoyancy production, I') call addfld('UW_sprod', ['ilev'], 'A', 'm2/s3', 'Shear production, I') call addfld('UW_kvh', ['ilev'], 'A', 'm2/s', 'Eddy diffusivity of heat, I') call addfld('UW_kvm', ['ilev'], 'A', 'm2/s', 'Eddy diffusivity of uv, I') call addfld('UW_pblh', horiz_only, 'A', 'm', 'PBLH, 1') call addfld('UW_pblhp', horiz_only, 'A', 'Pa', 'PBLH pressure, 1') call addfld('UW_tpert', horiz_only, 'A', 'K', 'Convective T excess, 1') call addfld('UW_qpert', horiz_only, 'A', 'kg/kg', 'Convective qt excess, I') call addfld('UW_wpert', horiz_only, 'A', 'm/s', 'Convective W excess, I') call addfld('UW_ustar', horiz_only, 'A', 'm/s', 'Surface Frictional Velocity, 1') call addfld('UW_tkes', horiz_only, 'A', 'm2/s2', 'Surface TKE, 1') call addfld('UW_minpblh',horiz_only, 'A', 'm', 'Minimum PBLH, 1') call addfld('UW_turbtype', ['ilev'], 'A', '1', 'Interface Turbulence Type, I') call addfld('UW_kbase_o', ['lev'], 'A', '1', 'Initial CL Base Exterbal Interface Index, CL') call addfld('UW_ktop_o', ['lev'], 'A', '1', 'Initial Top Exterbal Interface Index, CL') call addfld('UW_ncvfin_o',horiz_only,'A', '1', 'Initial Total Number of CL regimes, CL') call addfld('UW_kbase_mg', ['lev'], 'A', '1', 'kbase after merging, CL') call addfld('UW_ktop_mg', ['lev'], 'A', '1', 'ktop after merging, CL') call addfld('UW_ncvfin_mg',horiz_only,'A', '1', 'ncvfin after merging, CL') call addfld('UW_kbase_f', ['lev'], 'A', '1', 'Final kbase with SRCL, CL') call addfld('UW_ktop_f', ['lev'], 'A', '1', 'Final ktop with SRCL, CL') call addfld('UW_ncvfin_f',horiz_only,'A', '1', 'Final ncvfin with SRCL, CL') call addfld('UW_wet', ['lev'], 'A', 'm/s', 'Entrainment rate at CL top, CL') call addfld('UW_web', ['lev'], 'A', 'm/s', 'Entrainment rate at CL base, CL') call addfld('UW_jtbu', ['lev'], 'A', 'm/s2', 'Buoyancy jump across CL top, CL') call addfld('UW_jbbu', ['lev'], 'A', 'm/s2', 'Buoyancy jump across CL base, CL') call addfld('UW_evhc', ['lev'], 'A', '1', 'Evaporative enhancement factor, CL') call addfld('UW_jt2slv', ['lev'], 'A', 'J/kg', 'slv jump for evhc, CL') call addfld('UW_n2ht', ['lev'], 'A', 's-2', 'n2 at just below CL top interface, CL') call addfld('UW_n2hb', ['lev'], 'A', 's-2', 'n2 at just above CL base interface') call addfld('UW_lwp', ['lev'], 'A', 'kg/m2', 'LWP in the CL top layer, CL') call addfld('UW_optdepth', ['lev'], 'A', '1', 'Optical depth of the CL top layer, CL') call addfld('UW_radfrac', ['lev'], 'A', '1', 'Fraction of radiative cooling confined in the CL top') call addfld('UW_radf', ['lev'], 'A', 'm2/s3', 'Buoyancy production at the CL top by radf, I') call addfld('UW_wstar', ['lev'], 'A', 'm/s', 'Convective velocity, Wstar, CL') call addfld('UW_wstar3fact',['lev'], 'A', '1', 'Enhancement of wstar3 due to entrainment, CL') call addfld('UW_ebrk', ['lev'], 'A', 'm2/s2', 'CL-averaged TKE, CL') call addfld('UW_wbrk', ['lev'], 'A', 'm2/s2', 'CL-averaged W, CL') call addfld('UW_lbrk', ['lev'], 'A', 'm', 'CL internal thickness, CL') call addfld('UW_ricl', ['lev'], 'A', '1', 'CL-averaged Ri, CL') call addfld('UW_ghcl', ['lev'], 'A', '1', 'CL-averaged gh, CL') call addfld('UW_shcl', ['lev'], 'A', '1', 'CL-averaged sh, CL') call addfld('UW_smcl', ['lev'], 'A', '1', 'CL-averaged sm, CL') call addfld('UW_gh', ['ilev'], 'A', '1', 'gh at all interfaces, I') call addfld('UW_sh', ['ilev'], 'A', '1', 'sh at all interfaces, I') call addfld('UW_sm', ['ilev'], 'A', '1', 'sm at all interfaces, I') call addfld('UW_ria', ['ilev'], 'A', '1', 'ri at all interfaces, I') call addfld('UW_leng', ['ilev'], 'A', 'm/s', 'Turbulence length scale, I') ! For sedimentation-entrainment feedback analysis call addfld('UW_wsed', ['lev'], 'A', 'm/s', 'Sedimentation velocity at CL top, CL') end subroutine eddy_diff_init subroutine eddy_diff_tend(state, pbuf, cam_in, & ztodt, p, tint, rhoi, cldn, wstarent, & kvm_in, kvh_in, ksrftms, dragblj,tauresx, tauresy, & rrho, ustar, pblh, kvm, kvh, kvq, cgh, cgs, tpert, qpert, & tke, sprod, sfi, turbtype, sm_aw) use physics_types, only: physics_state use camsrfexch, only: cam_in_t use coords_1d, only: Coords1D type(physics_state), intent(in) :: state type(physics_buffer_desc), pointer, intent(in) :: pbuf(:) type(cam_in_t), intent(in) :: cam_in real(r8), intent(in) :: ztodt type(Coords1D), intent(in) :: p real(r8), intent(in) :: tint(pcols,pver+1) real(r8), intent(in) :: rhoi(pcols,pver+1) real(r8), intent(in) :: cldn(pcols,pver) logical, intent(in) :: wstarent real(r8), intent(in) :: kvm_in(pcols,pver+1) real(r8), intent(in) :: kvh_in(pcols,pver+1) real(r8), intent(in) :: ksrftms(pcols) real(r8), intent(in) :: dragblj(pcols,pver) ! Drag profile from Beljaars SGO form drag [ 1/s ] real(r8), intent(inout) :: tauresx(pcols) real(r8), intent(inout) :: tauresy(pcols) real(r8), intent(out) :: rrho(pcols) real(r8), intent(out) :: ustar(pcols) real(r8), intent(out) :: pblh(pcols) real(r8), intent(out) :: kvm(pcols,pver+1) real(r8), intent(out) :: kvh(pcols,pver+1) real(r8), intent(out) :: kvq(pcols,pver+1) real(r8), intent(out) :: cgh(pcols,pver+1) real(r8), intent(out) :: cgs(pcols,pver+1) real(r8), intent(out) :: tpert(pcols) real(r8), intent(out) :: qpert(pcols) real(r8), intent(out) :: tke(pcols,pver+1) real(r8), intent(out) :: sprod(pcols,pver+1) real(r8), intent(out) :: sfi(pcols,pver+1) integer(i4), intent(out) :: turbtype(pcols,pver+1) real(r8), intent(out) :: sm_aw(pcols,pver+1) integer :: i, k call compute_eddy_diff( pbuf, state%lchnk , & pcols , pver , state%ncol , state%t , tint, state%q(:,:,1) , ztodt , & state%q(:,:,ixcldliq) , state%q(:,:,ixcldice) , & state%s , p , rhoi, cldn , & state%zm , state%zi , state%pmid , state%pint , state%u , state%v , & cam_in%wsx, cam_in%wsy , cam_in%shf , cam_in%cflx(:,1) , wstarent , & rrho , ustar , pblh , kvm_in , kvh_in , kvm , & kvh , kvq , cgh , & cgs , tpert , qpert , tke , & sprod , sfi , & tauresx , tauresy , ksrftms , dragblj , turbtype , sm_aw ) ! The diffusivities from diag_TKE can be much larger than from HB in the free ! troposphere and upper atmosphere. These seem to be larger than observations, ! and in WACCM the gw_drag code is already applying an eddy diffusivity in the ! upper atmosphere. Optionally, adjust the diffusivities in the free troposphere ! or the upper atmosphere. ! ! NOTE: Further investigation should be done as to why the diffusivities are ! larger in diag_TKE. if ((kv_freetrop_scale /= 1._r8) .or. ((kv_top_scale /= 1._r8) .and. (kv_top_pressure > 0._r8))) then do i = 1, state%ncol do k = 1, pverp ! Outside of the boundary layer? if (state%zi(i,k) > pblh(i)) then ! In the upper atmosphere? if (state%pint(i,k) <= kv_top_pressure) then kvh(i,k) = kvh(i,k) * kv_top_scale kvm(i,k) = kvm(i,k) * kv_top_scale kvq(i,k) = kvq(i,k) * kv_top_scale else kvh(i,k) = kvh(i,k) * kv_freetrop_scale kvm(i,k) = kvm(i,k) * kv_freetrop_scale kvq(i,k) = kvq(i,k) * kv_freetrop_scale end if else exit end if end do end do end if end subroutine eddy_diff_tend !=============================================================================== ! ! ! !=============================================================================== ! subroutine compute_eddy_diff( pbuf, lchnk , & pcols , pver , ncol , t , tint, qv , ztodt , & ql , qi , s , p , rhoi, cldn , & z , zi , pmid , pi , u , v , & taux , tauy , shflx , qflx , wstarent , rrho , & ustar , pblh , kvm_in , kvh_in , kvm_out , kvh_out , kvq , & cgh , cgs , tpert , qpert , tke , & sprod , sfi , & tauresx, tauresy, ksrftms, dragblj, turbtype, sm_aw ) !-------------------------------------------------------------------- ! ! Purpose: Interface to compute eddy diffusivities. ! ! Eddy diffusivities are calculated in a fully implicit way ! ! through iteration process. ! ! Author: Sungsu Park. August. 2006. ! ! May. 2008. ! !-------------------------------------------------------------------- ! use diffusion_solver, only: compute_vdiff use cam_history, only: outfld use phys_debug_util, only: phys_debug_col use physconst, only: cpairv use pbl_utils, only: calc_ustar, austausch_atm use error_messages, only: handle_errmsg use coords_1d, only: Coords1D use wv_saturation, only: qsat use eddy_diff, only: trbintd, caleddy use physics_buffer, only: pbuf_get_field ! --------------- ! ! Input Variables ! ! --------------- ! type(physics_buffer_desc), pointer, intent(in) :: pbuf(:) integer, intent(in) :: lchnk integer, intent(in) :: pcols ! Number of atmospheric columns [ # ] integer, intent(in) :: pver ! Number of atmospheric layers [ # ] integer, intent(in) :: ncol ! Number of atmospheric columns [ # ] logical, intent(in) :: wstarent ! .true. means use the 'wstar' entrainment closure. real(r8), intent(in) :: ztodt ! Physics integration time step 2 delta-t [ s ] real(r8), intent(in) :: t(pcols,pver) ! Temperature [ K ] real(r8), intent(in) :: tint(pcols,pver+1) ! Temperature defined on interfaces [ K ] real(r8), intent(in) :: qv(pcols,pver) ! Water vapor specific humidity [ kg/kg ] real(r8), intent(in) :: ql(pcols,pver) ! Liquid water specific humidity [ kg/kg ] real(r8), intent(in) :: qi(pcols,pver) ! Ice specific humidity [ kg/kg ] real(r8), intent(in) :: s(pcols,pver) ! Dry static energy [ J/kg ] type(Coords1D), intent(in) :: p ! Pressure coordinates for solver [ Pa ] real(r8), intent(in) :: rhoi(pcols,pver+1) ! Density at interfaces [ kg/m3 ] real(r8), intent(in) :: cldn(pcols,pver) ! Stratiform cloud fraction [ fraction ] real(r8), intent(in) :: z(pcols,pver) ! Layer mid-point height above surface [ m ] real(r8), intent(in) :: zi(pcols,pver+1) ! Interface height above surface [ m ] real(r8), intent(in) :: pmid(pcols,pver) ! Layer mid-point pressure [ Pa ] real(r8), intent(in) :: pi(pcols,pver+1) ! Interface pressure [ Pa ] real(r8), intent(in) :: u(pcols,pver) ! Zonal velocity [ m/s ] real(r8), intent(in) :: v(pcols,pver) ! Meridional velocity [ m/s ] real(r8), intent(in) :: taux(pcols) ! Zonal wind stress at surface [ N/m2 ] real(r8), intent(in) :: tauy(pcols) ! Meridional wind stress at surface [ N/m2 ] real(r8), intent(in) :: shflx(pcols) ! Sensible heat flux at surface [ unit ? ] real(r8), intent(in) :: qflx(pcols) ! Water vapor flux at surface [ unit ? ] real(r8), intent(in) :: kvm_in(pcols,pver+1) ! kvm saved from last timestep [ m2/s ] real(r8), intent(in) :: kvh_in(pcols,pver+1) ! kvh saved from last timestep [ m2/s ] real(r8), intent(in) :: ksrftms(pcols) ! Surface drag coefficient of turbulent mountain stress [ unit ? ] real(r8), intent(in) :: dragblj(pcols,pver) ! Drag profile from Beljaars SGO form drag [ 1/s ] ! ---------------- ! ! Output Variables ! ! ---------------- ! real(r8), intent(out) :: kvm_out(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ] real(r8), intent(out) :: kvh_out(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ] real(r8), intent(out) :: kvq(pcols,pver+1) ! Eddy diffusivity for constituents, moisture and tracers [ m2/s ] ! (note not having '_out') real(r8), intent(out) :: rrho(pcols) ! Reciprocal of density at the lowest layer real(r8), intent(out) :: ustar(pcols) ! Surface friction velocity [ m/s ] real(r8), intent(out) :: pblh(pcols) ! PBL top height [ m ] real(r8), intent(out) :: cgh(pcols,pver+1) ! Counter-gradient term for heat [ J/kg/m ] real(r8), intent(out) :: cgs(pcols,pver+1) ! Counter-gradient star [ cg/flux ] real(r8), intent(out) :: tpert(pcols) ! Convective temperature excess [ K ] real(r8), intent(out) :: qpert(pcols) ! Convective humidity excess [ kg/kg ] real(r8), intent(out) :: tke(pcols,pver+1) ! Turbulent kinetic energy [ m2/s2 ] real(r8), intent(out) :: sprod(pcols,pver+1) ! Shear production [ m2/s3 ] real(r8), intent(out) :: sfi(pcols,pver+1) ! Interfacial layer saturation fraction [ fraction ] integer(i4), intent(out):: turbtype(pcols,pver+1) ! Turbulence type identifier at all interfaces [ no unit ] real(r8), intent(out) :: sm_aw(pcols,pver+1) ! Normalized Galperin instability function for momentum [ no unit ] ! This is 1 when neutral condition (Ri=0), ! 4.964 for maximum unstable case, and 0 when Ri > Ricrit=0.19. ! ---------------------- ! ! Input-Output Variables ! ! ---------------------- ! real(r8), intent(inout) :: tauresx(pcols) ! Residual stress to be added in vdiff to correct for turb real(r8), intent(inout) :: tauresy(pcols) ! Stress mismatch between sfc and atm accumulated in prior timesteps ! -------------- ! ! pbuf Variables ! ! -------------- ! real(r8), pointer :: qrl(:,:) ! LW radiative cooling rate real(r8), pointer :: wsedl(:,:) ! Sedimentation velocity ! of stratiform liquid cloud droplet [ m/s ] real(r8), pointer :: bprod(:,:) ! Buoyancy production of tke [ m2/s3 ] integer(i4), pointer :: ipbl(:) ! If 1, PBL is CL, while if 0, PBL is STL. integer(i4), pointer :: kpblh(:) ! Layer index containing PBL top within or at the base interface real(r8), pointer :: wstarPBL(:) ! Convective velocity within PBL [ m/s ] real(r8), pointer :: tkes(:) ! TKE at surface interface [ m2/s2 ] real(r8), pointer :: went(:) ! Entrainment rate at the PBL top interface [ m/s ] ! --------------- ! ! Local Variables ! ! --------------- ! integer icol integer i, k, iturb, status character(2048) :: warnstring ! Warning(s) to print character(128) :: errstring ! Error message real(r8) :: kvf(pcols,pver+1) ! Free atmospheric eddy diffusivity [ m2/s ] real(r8) :: kvm(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ] real(r8) :: kvh(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ] real(r8) :: kvm_preo(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ] real(r8) :: kvh_preo(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ] real(r8) :: kvm_pre(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ] real(r8) :: kvh_pre(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ] real(r8) :: errorPBL(pcols) ! Error function showing whether PBL produced convergent solution or not. ! [ unit ? ] real(r8) :: s2(pcols,pver) ! Shear squared, defined at interfaces except surface [ s-2 ] real(r8) :: n2(pcols,pver) ! Buoyancy frequency, defined at interfaces except surface [ s-2 ] real(r8) :: ri(pcols,pver) ! Richardson number, 'n2/s2', defined at interfaces except surface [ s-2 ] real(r8) :: pblhp(pcols) ! PBL top pressure [ Pa ] real(r8) :: minpblh(pcols) ! Minimum PBL height based on surface stress real(r8) :: qt(pcols,pver) ! Total specific humidity [ kg/kg ] real(r8) :: sfuh(pcols,pver) ! Saturation fraction in upper half-layer [ fraction ] real(r8) :: sflh(pcols,pver) ! Saturation fraction in lower half-layer [ fraction ] real(r8) :: sl(pcols,pver) ! Liquid water static energy [ J/kg ] real(r8) :: slv(pcols,pver) ! Liquid water virtual static energy [ J/kg ] real(r8) :: slslope(pcols,pver) ! Slope of 'sl' in each layer real(r8) :: qtslope(pcols,pver) ! Slope of 'qt' in each layer real(r8) :: qvfd(pcols,pver) ! Specific humidity for diffusion [ kg/kg ] real(r8) :: tfd(pcols,pver) ! Temperature for diffusion [ K ] real(r8) :: slfd(pcols,pver) ! Liquid static energy [ J/kg ] real(r8) :: qtfd(pcols,pver) ! Total specific humidity [ kg/kg ] real(r8) :: qlfd(pcols,pver) ! Liquid water specific humidity for diffusion [ kg/kg ] real(r8) :: ufd(pcols,pver) ! U-wind for diffusion [ m/s ] real(r8) :: vfd(pcols,pver) ! V-wind for diffusion [ m/s ] ! Buoyancy coefficients : w'b' = ch * w'sl' + cm * w'qt' real(r8) :: chu(pcols,pver+1) ! Heat buoyancy coef for dry states, defined at each interface, finally. real(r8) :: chs(pcols,pver+1) ! Heat buoyancy coef for sat states, defined at each interface, finally. real(r8) :: cmu(pcols,pver+1) ! Moisture buoyancy coef for dry states, ! defined at each interface, finally. real(r8) :: cms(pcols,pver+1) ! Moisture buoyancy coef for sat states, ! defined at each interface, finally. real(r8) :: jnk1d(pcols) real(r8) :: jnk2d(pcols,pver+1) real(r8) :: zero(pcols) real(r8) :: zero2d(pcols,pver+1) real(r8) :: es ! Saturation vapor pressure real(r8) :: qs ! Saturation specific humidity real(r8) :: ep2, templ, temps ! ------------------------------- ! ! Variables for diagnostic output ! ! ------------------------------- ! real(r8) :: wpert(pcols) ! Turbulent velocity excess [ m/s ] real(r8) :: kbase_o(pcols,ncvmax) ! Original external base interface index of CL from 'exacol' real(r8) :: ktop_o(pcols,ncvmax) ! Original external top interface index of CL from 'exacol' real(r8) :: ncvfin_o(pcols) ! Original number of CLs from 'exacol' real(r8) :: kbase_mg(pcols,ncvmax) ! 'kbase' after extending-merging from 'zisocl' real(r8) :: ktop_mg(pcols,ncvmax) ! 'ktop' after extending-merging from 'zisocl' real(r8) :: ncvfin_mg(pcols) ! 'ncvfin' after extending-merging from 'zisocl' real(r8) :: kbase_f(pcols,ncvmax) ! Final 'kbase' after extending-merging & including SRCL real(r8) :: ktop_f(pcols,ncvmax) ! Final 'ktop' after extending-merging & including SRCL real(r8) :: ncvfin_f(pcols) ! Final 'ncvfin' after extending-merging & including SRCL real(r8) :: wet(pcols,ncvmax) ! Entrainment rate at the CL top [ m/s ] real(r8) :: web(pcols,ncvmax) ! Entrainment rate at the CL base [ m/s ]. ! Set to zero if CL is based at surface. real(r8) :: jtbu(pcols,ncvmax) ! Buoyancy jump across the CL top [ m/s2 ] real(r8) :: jbbu(pcols,ncvmax) ! Buoyancy jump across the CL base [ m/s2 ] real(r8) :: evhc(pcols,ncvmax) ! Evaporative enhancement factor at the CL top real(r8) :: jt2slv(pcols,ncvmax) ! Jump of slv ( across two layers ) at CL top used only for evhc [ J/kg ] real(r8) :: n2ht(pcols,ncvmax) ! n2 defined at the CL top interface but using ! sfuh(kt) instead of sfi(kt) [ s-2 ] real(r8) :: n2hb(pcols,ncvmax) ! n2 defined at the CL base interface but using ! sflh(kb-1) instead of sfi(kb) [ s-2 ] real(r8) :: lwp(pcols,ncvmax) ! LWP in the CL top layer [ kg/m2 ] real(r8) :: opt_depth(pcols,ncvmax) ! Optical depth of the CL top layer real(r8) :: radinvfrac(pcols,ncvmax) ! Fraction of radiative cooling confined in the top portion of CL top layer real(r8) :: radf(pcols,ncvmax) ! Buoyancy production at the CL top due to LW radiative cooling [ m2/s3 ] real(r8) :: wstar(pcols,ncvmax) ! Convective velocity in each CL [ m/s ] real(r8) :: wstar3fact(pcols,ncvmax) ! Enhancement of 'wstar3' due to entrainment (inverse) [ no unit ] real(r8) :: ebrk(pcols,ncvmax) ! Net mean TKE of CL including entrainment effect [ m2/s2 ] real(r8) :: wbrk(pcols,ncvmax) ! Net mean normalized TKE (W) of CL, ! 'ebrk/b1' including entrainment effect [ m2/s2 ] real(r8) :: lbrk(pcols,ncvmax) ! Energetic internal thickness of CL [m] real(r8) :: ricl(pcols,ncvmax) ! CL internal mean Richardson number real(r8) :: ghcl(pcols,ncvmax) ! Half of normalized buoyancy production of CL real(r8) :: shcl(pcols,ncvmax) ! Galperin instability function of heat-moisture of CL real(r8) :: smcl(pcols,ncvmax) ! Galperin instability function of mementum of CL real(r8) :: ghi(pcols,pver+1) ! Half of normalized buoyancy production at all interfaces real(r8) :: shi(pcols,pver+1) ! Galperin instability function of heat-moisture at all interfaces real(r8) :: smi(pcols,pver+1) ! Galperin instability function of heat-moisture at all interfaces real(r8) :: rii(pcols,pver+1) ! Interfacial Richardson number defined at all interfaces real(r8) :: lengi(pcols,pver+1) ! Turbulence length scale at all interfaces [ m ] real(r8) :: wcap(pcols,pver+1) ! Normalized TKE at all interfaces [ m2/s2 ] ! For sedimentation-entrainment feedback real(r8) :: wsed(pcols,ncvmax) ! Sedimentation velocity at the top of each CL [ m/s ] ! ---------- ! ! Parameters ! ! ---------- ! logical, parameter :: use_kvf = .false. ! .true. (.false.) : initialize kvh/kvm = kvf ( 0. ) real(r8), parameter :: lambda = 0.5_r8 ! Under-relaxation factor ( 0 < lambda =< 1 ) ! ---------- ! ! Initialize ! ! ---------- ! zero(:) = 0._r8 zero2d(:,:) = 0._r8 ! ---------------------------------------------- ! ! Get LW radiative heating out of physics buffer ! ! ---------------------------------------------- ! call pbuf_get_field(pbuf, qrl_idx, qrl) call pbuf_get_field(pbuf, wsedl_idx, wsedl) ! These fields are put into the pbuf for UNICON only. if (unicon_is_on) then call pbuf_get_field(pbuf, bprod_idx, bprod) call pbuf_get_field(pbuf, ipbl_idx, ipbl) call pbuf_get_field(pbuf, kpblh_idx, kpblh) call pbuf_get_field(pbuf, wstarPBL_idx, wstarPBL) call pbuf_get_field(pbuf, tkes_idx, tkes) call pbuf_get_field(pbuf, went_idx, went) else allocate(bprod(pcols,pverp), ipbl(pcols), kpblh(pcols), wstarPBL(pcols), tkes(pcols), went(pcols)) end if ! ----------------------- ! ! Main Computation Begins ! ! ----------------------- ! ufd(:ncol,:) = u(:ncol,:) vfd(:ncol,:) = v(:ncol,:) tfd(:ncol,:) = t(:ncol,:) qvfd(:ncol,:) = qv(:ncol,:) qlfd(:ncol,:) = ql(:ncol,:) 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 ) ! 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 == 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 if (use_kvf) then call austausch_atm(pcols, ncol, pver, ntop_eddy, nbot_eddy, & ml2, ri, s2, kvf ) else kvf = 0._r8 end if ! Initialize kvh/kvm to send to caleddy, depending on model timestep and iteration number ! This is necessary for 'wstar-based' entrainment closure. if( iturb == 1 ) then if( is_first_step() ) then ! First iteration of first model timestep : Use free tropospheric value or zero. kvh(:ncol,:) = kvf(:ncol,:) kvm(:ncol,:) = kvf(:ncol,:) 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 , went , 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 , wsed, & warnstring, errstring) if (trim(warnstring) /= "") then write(iulog,*) "eddy_diff_cam: Messages from caleddy follow." write(iulog,*) warnstring end if call handle_errmsg(errstring, subname="caleddy") ! Calculate errorPBL to check whether PBL produced convergent solutions or not. if( iturb == 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 > 1 .and. iturb < 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 < 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,:) ! 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 , tint, & p , t , rhoi, ztodt , taux , & tauy , shflx , qflx , & kvh_out , kvm_out , kvh_out , cgs , cgh , & zi , ksrftms , dragblj , & zero , fieldlist_wet, fieldlist_molec, & ufd , vfd , qtfd , slfd , & jnk1d , jnk1d , jnk2d , jnk1d , errstring , & tauresx , tauresy , 0 , cpairv(:,:,lchnk), zero, & .false., .false. ) call handle_errmsg(errstring, subname="compute_vdiff", & extra_msg="compute_vdiff called from eddy_diff_cam") ! 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 be 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) - gravit*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) + (latvap+latice) * qi(i,k) - gravit*z(i,k)) / cpair end do end do endif end do ! End of 'iturb' iteration kvq(:ncol,:) = kvh_out(:ncol,:) ! Compute 'wstar' within the PBL for use in the future convection scheme. do i = 1, ncol if(ipbl(i) == 1) 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( 'WGUSTD' , wpert, pcols, lchnk ) call outfld( 'BPROD ', bprod, pcols, lchnk ) call outfld( 'SPROD ', sprod, pcols, lchnk ) call outfld( 'SFI ', sfi, pcols, lchnk ) 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 ) call outfld( 'UW_wsed', wsed, pcols, lchnk ) if (.not. unicon_is_on) then deallocate(bprod, ipbl, kpblh, wstarPBL, tkes, went) end if end subroutine compute_eddy_diff end module eddy_diff_cam