LCOV - code coverage report
Current view: | top level - physics/cam - eddy_diff_cam.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 393 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 5 | 0.0 %

              Line data    Source code
           1             : module eddy_diff_cam
           2             :
           3             : use shr_kind_mod, only: i4 => shr_kind_i4, r8 => shr_kind_r8
           4             : use ppgrid, only: pcols, pver, pverp
           5             : use cam_logfile, only: iulog
           6             : use cam_abortutils, only: endrun
           7             : use physconst, only: gravit, cpair, rair, zvir, latvap, latice, karman
           8             : use diffusion_solver, only: vdiff_selector
           9             : use eddy_diff, only: ncvmax
          10             : use time_manager, only: is_first_step
          11             : use physics_buffer, only: physics_buffer_desc
          12             : use spmd_utils, only: masterproc
          13             : use phys_control, only: phys_getopts
          14             :
          15             : implicit none
          16             : private
          17             :
          18             : public :: eddy_diff_readnl
          19             : public :: eddy_diff_register
          20             : public :: eddy_diff_init
          21             : public :: eddy_diff_tend
          22             :
          23             : ! Is UNICON switched on (and thus interacting with eddy_diff via pbuf)?
          24             : logical :: unicon_is_on
          25             :
          26             : ! Number of iterations for solution
          27             : integer, parameter :: nturb = 5
          28             :
          29             : ! Logical switches for moist mixing ratio diffusion
          30             : type(vdiff_selector) :: fieldlist_wet
          31             : ! Logical switches for molecular diffusion
          32             : ! (Molecular diffusion is not done here.)
          33             : type(vdiff_selector) :: fieldlist_molec
          34             :
          35             : integer :: ntop_eddy, nbot_eddy
          36             :
          37             : ! Cloud mass constituent indices
          38             : integer :: ixcldliq, ixcldice
          39             :
          40             : ! input pbuf field indices
          41             : integer :: qrl_idx   = -1
          42             : integer :: wsedl_idx = -1
          43             :
          44             : ! output pbuf field indices for UNICON
          45             : integer :: bprod_idx    = -1
          46             : integer :: ipbl_idx     = -1
          47             : integer :: kpblh_idx    = -1
          48             : integer :: wstarPBL_idx = -1
          49             : integer :: tkes_idx     = -1
          50             : integer :: went_idx     = -1
          51             :
          52             : ! Mixing lengths squared.
          53             : ! Used for computing free air diffusivity.
          54             : real(r8) :: ml2(pver+1)
          55             :
          56             : ! Various namelist options to limit or tweak the effects of eddy diffusion.
          57             :
          58             : ! Pressure defining the bottom of the upper atmosphere for kvh scaling (Pa)
          59             : real(r8) :: kv_top_pressure = 0._r8
          60             : ! Eddy diffusivity scale factor for upper atmosphere
          61             : real(r8) :: kv_top_scale = 1._r8
          62             : ! Eddy diffusivity scale factor for the free troposphere
          63             : real(r8) :: kv_freetrop_scale = 1._r8
          64             :
          65             : ! The following all have to be set in all cases.
          66             : real(r8), parameter  :: unset_r8 = huge(1._r8)
          67             : ! Maximum master length for diag_TKE
          68             : real(r8) :: eddy_lbulk_max = unset_r8
          69             : ! Maximum dissipation length for diag_TKE
          70             : real(r8) :: eddy_leng_max = unset_r8
          71             : ! Bottom pressure level (hPa) for eddy_leng_max
          72             : real(r8) :: eddy_max_bot_pressure = unset_r8
          73             : ! Moist entrainment enhancement param
          74             : real(r8) :: eddy_moist_entrain_a2l = unset_r8
          75             :
          76             : contains
          77             :
          78           0 : subroutine eddy_diff_readnl(nlfile)
          79             :   use namelist_utils, only: find_group_name
          80             :   use units, only: getunit, freeunit
          81             :   use spmd_utils, only: masterprocid, mpi_real8, mpicom
          82             :   use shr_log_mod, only: errMsg => shr_log_errMsg
          83             :
          84             :   ! filepath for file containing namelist input
          85             :   character(len=*), intent(in) :: nlfile
          86             :
          87             :   ! file unit and error code
          88             :   integer :: unitn, ierr
          89             :
          90             :   character(len=*), parameter :: subname = 'eddy_diff_readnl'
          91             :
          92             :   namelist /eddy_diff_nl/ kv_top_pressure, kv_top_scale, &
          93             :        kv_freetrop_scale, eddy_lbulk_max, eddy_leng_max, &
          94             :        eddy_max_bot_pressure, eddy_moist_entrain_a2l
          95             :
          96           0 :   if (masterproc) then
          97           0 :      unitn = getunit()
          98           0 :      open( unitn, file=trim(nlfile), status='old' )
          99           0 :      call find_group_name(unitn, 'eddy_diff_nl', status=ierr)
         100           0 :      if (ierr == 0) then
         101           0 :         read(unitn, eddy_diff_nl, iostat=ierr)
         102             :      end if
         103           0 :      if (ierr /= 0) then
         104           0 :         call endrun(subname // ':: ERROR reading namelist')
         105             :      end if
         106           0 :      close(unitn)
         107           0 :      call freeunit(unitn)
         108             :   end if
         109             :
         110           0 :   call mpi_bcast(kv_top_pressure,   1,   mpi_real8, masterprocid, mpicom, ierr)
         111           0 :   if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
         112           0 :   call mpi_bcast(kv_top_scale,      1,   mpi_real8, masterprocid, mpicom, ierr)
         113           0 :   if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
         114           0 :   call mpi_bcast(kv_freetrop_scale, 1,   mpi_real8, masterprocid, mpicom, ierr)
         115           0 :   if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
         116             :
         117           0 :   call mpi_bcast(eddy_lbulk_max,    1,   mpi_real8, masterprocid, mpicom, ierr)
         118           0 :   if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
         119           0 :   call mpi_bcast(eddy_leng_max,     1,   mpi_real8, masterprocid, mpicom, ierr)
         120           0 :   if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
         121           0 :   call mpi_bcast(eddy_max_bot_pressure,     1,   mpi_real8, masterprocid, mpicom, ierr)
         122           0 :   if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
         123           0 :   call mpi_bcast(eddy_moist_entrain_a2l,    1,   mpi_real8, masterprocid, mpicom, ierr)
         124           0 :   if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
         125             :
         126           0 : end subroutine eddy_diff_readnl
         127             :
         128           0 : subroutine eddy_diff_register()
         129             :   use physics_buffer, only: pbuf_add_field, dtype_r8, dtype_i4
         130             :
         131             :   character(len=16) :: shallow_scheme
         132             :
         133             :   ! Check for UNICON and add relevant pbuf entries.
         134           0 :   call phys_getopts(shallow_scheme_out=shallow_scheme)
         135             :
         136           0 :   unicon_is_on = (shallow_scheme == "UNICON")
         137             :
         138           0 :   if (unicon_is_on) then
         139           0 :      call pbuf_add_field('bprod',    'global', dtype_r8, (/pcols,pverp/), bprod_idx)
         140           0 :      call pbuf_add_field('ipbl',     'global', dtype_i4, (/pcols/),       ipbl_idx)
         141           0 :      call pbuf_add_field('kpblh',    'global', dtype_i4, (/pcols/),       kpblh_idx)
         142           0 :      call pbuf_add_field('wstarPBL', 'global', dtype_r8, (/pcols/),       wstarPBL_idx)
         143           0 :      call pbuf_add_field('tkes',     'global', dtype_r8, (/pcols/),       tkes_idx)
         144           0 :      call pbuf_add_field('went',     'global', dtype_r8, (/pcols/),       went_idx)
         145             :   end if
         146             :
         147           0 : end subroutine eddy_diff_register
         148             :
         149           0 : subroutine eddy_diff_init(pbuf2d, ntop_eddy_in, nbot_eddy_in)
         150             :
         151           0 :   use error_messages, only: handle_errmsg
         152             :   use cam_history, only: addfld, add_default, horiz_only
         153             :   use constituents, only: cnst_get_ind
         154             :   use ref_pres, only: pref_mid
         155             :   use diffusion_solver, only: new_fieldlist_vdiff, vdiff_select
         156             :   use eddy_diff, only: init_eddy_diff
         157             :   use physics_buffer, only: pbuf_set_field, pbuf_get_index
         158             :
         159             :   type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! Physics buffer
         160             :   integer,  intent(in) :: ntop_eddy_in ! Top interface level to which eddy vertical diffusivity is applied ( = 1 )
         161             :   integer,  intent(in) :: nbot_eddy_in ! Bottom interface level to which eddy vertical diffusivity is applied ( = pver )
         162             :
         163             :   character(len=128) :: errstring
         164             :
         165             :   real(r8) :: leng_max(pver)
         166             :   integer :: k
         167             :
         168             :   logical :: history_amwg
         169             :
         170           0 :   ntop_eddy = ntop_eddy_in
         171           0 :   nbot_eddy = nbot_eddy_in
         172             :
         173           0 :   do k = 1, pver
         174           0 :      if (pref_mid(k) <= eddy_max_bot_pressure*1.e2_r8) then
         175           0 :         leng_max(k) = eddy_leng_max
         176             :      else
         177           0 :         leng_max(k) = 40.e3_r8
         178             :      end if
         179             :   end do
         180             :
         181           0 :   if (masterproc) then
         182           0 :      write(iulog,*)'init_eddy_diff: nturb=',nturb
         183           0 :      write(iulog,*)'init_eddy_diff: eddy_leng_max=',eddy_leng_max,' lbulk_max=',eddy_lbulk_max
         184           0 :      do k = 1,pver
         185           0 :         write(iulog,*)'init_eddy_diff:',k,pref_mid(k),'leng_max=',leng_max(k)
         186             :      end do
         187             :   end if
         188             :
         189             :   call init_eddy_diff(pver, gravit, cpair, rair, zvir, &
         190             :        latvap, latice, ntop_eddy, nbot_eddy, karman, &
         191             :        eddy_lbulk_max, leng_max, &
         192           0 :        eddy_moist_entrain_a2l, errstring)
         193             :
         194           0 :   call handle_errmsg(errstring, subname="init_eddy_diff")
         195             :
         196             :   ! Set the square of the mixing lengths.
         197           0 :   ml2(1:ntop_eddy) = 0._r8
         198           0 :   do k = ntop_eddy + 1, nbot_eddy
         199           0 :      ml2(k) = 30.0_r8**2
         200             :   end do
         201           0 :   ml2(nbot_eddy+1:pver+1) = 0._r8
         202             :
         203             :   ! Get fieldlists to pass to diffusion solver.
         204           0 :   fieldlist_wet   = new_fieldlist_vdiff(1)
         205           0 :   fieldlist_molec = new_fieldlist_vdiff(1)
         206             :
         207             :   call handle_errmsg(vdiff_select(fieldlist_wet,'s'), &
         208           0 :        subname="vdiff_select")
         209             :   call handle_errmsg(vdiff_select(fieldlist_wet,'q',1), &
         210           0 :        subname="vdiff_select")
         211             :   call handle_errmsg(vdiff_select(fieldlist_wet,'u'), &
         212           0 :        subname="vdiff_select")
         213             :   call handle_errmsg(vdiff_select(fieldlist_wet,'v'), &
         214           0 :        subname="vdiff_select")
         215             :
         216             :   ! Cloud mass constituents
         217           0 :   call cnst_get_ind('CLDLIQ', ixcldliq)
         218           0 :   call cnst_get_ind('CLDICE', ixcldice)
         219             :
         220             :   ! Input pbuf fields
         221           0 :   qrl_idx   = pbuf_get_index('QRL')
         222           0 :   wsedl_idx = pbuf_get_index('WSEDL')
         223             :
         224             :   ! Initialize output pbuf fields
         225           0 :   if (is_first_step() .and. unicon_is_on) then
         226           0 :      call pbuf_set_field(pbuf2d, bprod_idx,    1.0e-5_r8)
         227           0 :      call pbuf_set_field(pbuf2d, ipbl_idx,     0    )
         228           0 :      call pbuf_set_field(pbuf2d, kpblh_idx,    1    )
         229           0 :      call pbuf_set_field(pbuf2d, wstarPBL_idx, 0.0_r8)
         230           0 :      call pbuf_set_field(pbuf2d, tkes_idx,     0.0_r8)
         231           0 :      call pbuf_set_field(pbuf2d, went_idx,     0.0_r8)
         232             :   end if
         233             :
         234             :   ! Scheme-specific default output.
         235           0 :   call phys_getopts(history_amwg_out=history_amwg)
         236             :
         237           0 :   call addfld('WGUSTD', horiz_only,  'A',          'm/s',  'wind gusts from turbulence'                            )
         238           0 :   if (history_amwg) then
         239           0 :      call add_default( 'WGUSTD  ', 1, ' ' )
         240             :   end if
         241             :
         242             :   ! ------------------------------------------------------------------- !
         243             :   ! Writing outputs for detailed analysis of UW moist turbulence scheme !
         244             :   ! ------------------------------------------------------------------- !
         245             :
         246           0 :   call addfld( 'BPROD',   ['ilev'],  'A',        'm2/s3',  'Buoyancy Production'                                   )
         247           0 :   call addfld( 'SFI',     ['ilev'],  'A',            '1',  'Interface-layer sat frac'                              )
         248           0 :   call addfld( 'SPROD',   ['ilev'],  'A',        'm2/s3',  'Shear Production'                                      )
         249             :
         250             :
         251           0 :   call addfld('UW_errorPBL',horiz_only,'A',         'm2/s',  'Error function of UW PBL')
         252           0 :   call addfld('UW_n2',        ['lev'], 'A',          's-2',  'Buoyancy Frequency, LI')
         253           0 :   call addfld('UW_s2',        ['lev'], 'A',          's-2',  'Shear Frequency, LI')
         254           0 :   call addfld('UW_ri',        ['lev'], 'A',            '1',  'Interface Richardson Number, I')
         255           0 :   call addfld('UW_sfuh',      ['lev'], 'A',            '1',  'Upper-Half Saturation Fraction, L')
         256           0 :   call addfld('UW_sflh',      ['lev'], 'A',            '1',  'Lower-Half Saturation Fraction, L')
         257           0 :   call addfld('UW_sfi',      ['ilev'], 'A',            '1',  'Interface Saturation Fraction, I')
         258           0 :   call addfld('UW_cldn',      ['lev'], 'A',            '1',  'Cloud Fraction, L')
         259           0 :   call addfld('UW_qrl',       ['lev'], 'A', 'gravity W/m2',  'LW cooling rate, L')
         260           0 :   call addfld('UW_ql',        ['lev'], 'A',        'kg/kg',  'ql(LWC), L')
         261           0 :   call addfld('UW_chu',      ['ilev'], 'A', 'gravity kg/J',  'Buoyancy Coefficient, chu, I')
         262           0 :   call addfld('UW_chs',      ['ilev'], 'A', 'gravity kg/J',  'Buoyancy Coefficient, chs, I')
         263           0 :   call addfld('UW_cmu',      ['ilev'], 'A','gravity/kg/kg',  'Buoyancy Coefficient, cmu, I')
         264           0 :   call addfld('UW_cms',      ['ilev'], 'A','gravity/kg/kg',  'Buoyancy Coefficient, cms, I')
         265           0 :   call addfld('UW_tke',      ['ilev'], 'A',        'm2/s2',  'TKE, I')
         266           0 :   call addfld('UW_wcap',     ['ilev'], 'A',        'm2/s2',  'Wcap, I')
         267           0 :   call addfld('UW_bprod',    ['ilev'], 'A',        'm2/s3',  'Buoyancy production, I')
         268           0 :   call addfld('UW_sprod',    ['ilev'], 'A',        'm2/s3',  'Shear production, I')
         269           0 :   call addfld('UW_kvh',      ['ilev'], 'A',         'm2/s',  'Eddy diffusivity of heat, I')
         270           0 :   call addfld('UW_kvm',      ['ilev'], 'A',         'm2/s',  'Eddy diffusivity of uv, I')
         271           0 :   call addfld('UW_pblh',   horiz_only, 'A',            'm',  'PBLH, 1')
         272           0 :   call addfld('UW_pblhp',  horiz_only, 'A',           'Pa',  'PBLH pressure, 1')
         273           0 :   call addfld('UW_tpert',  horiz_only, 'A',            'K',  'Convective T excess, 1')
         274           0 :   call addfld('UW_qpert',  horiz_only, 'A',        'kg/kg',  'Convective qt excess, I')
         275           0 :   call addfld('UW_wpert',  horiz_only, 'A',          'm/s',  'Convective W excess, I')
         276           0 :   call addfld('UW_ustar',  horiz_only, 'A',          'm/s',  'Surface Frictional Velocity, 1')
         277           0 :   call addfld('UW_tkes',   horiz_only, 'A',        'm2/s2',  'Surface TKE, 1')
         278           0 :   call addfld('UW_minpblh',horiz_only, 'A',            'm',  'Minimum PBLH, 1')
         279           0 :   call addfld('UW_turbtype', ['ilev'], 'A',            '1',  'Interface Turbulence Type, I')
         280           0 :   call addfld('UW_kbase_o',   ['lev'], 'A',            '1',  'Initial CL Base Exterbal Interface Index, CL')
         281           0 :   call addfld('UW_ktop_o',    ['lev'], 'A',            '1',  'Initial Top Exterbal Interface Index, CL')
         282           0 :   call addfld('UW_ncvfin_o',horiz_only,'A',            '1',  'Initial Total Number of CL regimes, CL')
         283           0 :   call addfld('UW_kbase_mg',  ['lev'], 'A',            '1',  'kbase after merging, CL')
         284           0 :   call addfld('UW_ktop_mg',   ['lev'], 'A',            '1',  'ktop after merging, CL')
         285           0 :   call addfld('UW_ncvfin_mg',horiz_only,'A',           '1',  'ncvfin after merging, CL')
         286           0 :   call addfld('UW_kbase_f',   ['lev'], 'A',            '1',  'Final kbase with SRCL, CL')
         287           0 :   call addfld('UW_ktop_f',    ['lev'], 'A',            '1',  'Final ktop with SRCL, CL')
         288           0 :   call addfld('UW_ncvfin_f',horiz_only,'A',            '1',  'Final ncvfin with SRCL, CL')
         289           0 :   call addfld('UW_wet',       ['lev'], 'A',          'm/s',  'Entrainment rate at CL top, CL')
         290           0 :   call addfld('UW_web',       ['lev'], 'A',          'm/s',  'Entrainment rate at CL base, CL')
         291           0 :   call addfld('UW_jtbu',      ['lev'], 'A',         'm/s2',  'Buoyancy jump across CL top, CL')
         292           0 :   call addfld('UW_jbbu',      ['lev'], 'A',         'm/s2',  'Buoyancy jump across CL base, CL')
         293           0 :   call addfld('UW_evhc',      ['lev'], 'A',            '1',  'Evaporative enhancement factor, CL')
         294           0 :   call addfld('UW_jt2slv',    ['lev'], 'A',         'J/kg',  'slv jump for evhc, CL')
         295           0 :   call addfld('UW_n2ht',      ['lev'], 'A',          's-2',  'n2 at just below CL top interface, CL')
         296           0 :   call addfld('UW_n2hb',      ['lev'], 'A',          's-2',  'n2 at just above CL base interface')
         297           0 :   call addfld('UW_lwp',       ['lev'], 'A',        'kg/m2',  'LWP in the CL top layer, CL')
         298           0 :   call addfld('UW_optdepth',  ['lev'], 'A',            '1',  'Optical depth of the CL top layer, CL')
         299           0 :   call addfld('UW_radfrac',   ['lev'], 'A',            '1',  'Fraction of radiative cooling confined in the CL top')
         300           0 :   call addfld('UW_radf',      ['lev'], 'A',        'm2/s3',  'Buoyancy production at the CL top by radf, I')
         301           0 :   call addfld('UW_wstar',     ['lev'], 'A',          'm/s',  'Convective velocity, Wstar, CL')
         302           0 :   call addfld('UW_wstar3fact',['lev'], 'A',            '1',  'Enhancement of wstar3 due to entrainment, CL')
         303           0 :   call addfld('UW_ebrk',      ['lev'], 'A',        'm2/s2',  'CL-averaged TKE, CL')
         304           0 :   call addfld('UW_wbrk',      ['lev'], 'A',        'm2/s2',  'CL-averaged W, CL')
         305           0 :   call addfld('UW_lbrk',      ['lev'], 'A',            'm',  'CL internal thickness, CL')
         306           0 :   call addfld('UW_ricl',      ['lev'], 'A',            '1',  'CL-averaged Ri, CL')
         307           0 :   call addfld('UW_ghcl',      ['lev'], 'A',            '1',  'CL-averaged gh, CL')
         308           0 :   call addfld('UW_shcl',      ['lev'], 'A',            '1',  'CL-averaged sh, CL')
         309           0 :   call addfld('UW_smcl',      ['lev'], 'A',            '1',  'CL-averaged sm, CL')
         310           0 :   call addfld('UW_gh',       ['ilev'], 'A',            '1',  'gh at all interfaces, I')
         311           0 :   call addfld('UW_sh',       ['ilev'], 'A',            '1',  'sh at all interfaces, I')
         312           0 :   call addfld('UW_sm',       ['ilev'], 'A',            '1',  'sm at all interfaces, I')
         313           0 :   call addfld('UW_ria',      ['ilev'], 'A',            '1',  'ri at all interfaces, I')
         314           0 :   call addfld('UW_leng',     ['ilev'], 'A',          'm/s',  'Turbulence length scale, I')
         315             :   ! For sedimentation-entrainment feedback analysis
         316           0 :   call addfld('UW_wsed',      ['lev'], 'A',          'm/s',  'Sedimentation velocity at CL top, CL')
         317             :
         318           0 : end subroutine eddy_diff_init
         319             :
         320           0 : subroutine eddy_diff_tend(state, pbuf, cam_in, &
         321             :      ztodt, p, tint, rhoi, cldn, wstarent, &
         322             :      kvm_in, kvh_in, ksrftms, dragblj,tauresx, tauresy, &
         323             :      rrho, ustar, pblh, kvm, kvh, kvq, cgh, cgs, tpert, qpert, &
         324             :      tke, sprod, sfi)
         325             :
         326           0 :   use physics_types, only: physics_state
         327             :   use camsrfexch, only: cam_in_t
         328             :   use coords_1d, only: Coords1D
         329             :
         330             :   type(physics_state), intent(in) :: state
         331             :   type(physics_buffer_desc), pointer, intent(in) :: pbuf(:)
         332             :   type(cam_in_t), intent(in) :: cam_in
         333             :   real(r8), intent(in) :: ztodt
         334             :   type(Coords1D), intent(in) :: p
         335             :   real(r8), intent(in) :: tint(pcols,pver+1)
         336             :   real(r8), intent(in) :: rhoi(pcols,pver+1)
         337             :   real(r8), intent(in) :: cldn(pcols,pver)
         338             :   logical, intent(in) :: wstarent
         339             :   real(r8), intent(in) :: kvm_in(pcols,pver+1)
         340             :   real(r8), intent(in) :: kvh_in(pcols,pver+1)
         341             :   real(r8), intent(in) :: ksrftms(pcols)
         342             :   real(r8), intent(in) :: dragblj(pcols,pver)       ! Drag profile from Beljaars SGO form drag [ 1/s ]
         343             :   real(r8), intent(inout) :: tauresx(pcols)
         344             :   real(r8), intent(inout) :: tauresy(pcols)
         345             :   real(r8), intent(out) :: rrho(pcols)
         346             :   real(r8), intent(out) :: ustar(pcols)
         347             :   real(r8), intent(out) :: pblh(pcols)
         348             :   real(r8), intent(out) :: kvm(pcols,pver+1)
         349             :   real(r8), intent(out) :: kvh(pcols,pver+1)
         350             :   real(r8), intent(out) :: kvq(pcols,pver+1)
         351             :   real(r8), intent(out) :: cgh(pcols,pver+1)
         352             :   real(r8), intent(out) :: cgs(pcols,pver+1)
         353             :   real(r8), intent(out) :: tpert(pcols)
         354             :   real(r8), intent(out) :: qpert(pcols)
         355             :   real(r8), intent(out) :: tke(pcols,pver+1)
         356             :   real(r8), intent(out) :: sprod(pcols,pver+1)
         357             :   real(r8), intent(out) :: sfi(pcols,pver+1)
         358             :
         359             :   integer :: i, k
         360             :
         361             :   call compute_eddy_diff( pbuf, state%lchnk    ,                                     &
         362             :        pcols    , pver        , state%ncol       , state%t    , tint, state%q(:,:,1) , ztodt   , &
         363             :        state%q(:,:,ixcldliq)  , state%q(:,:,ixcldice)   ,                            &
         364             :        state%s  , p           , rhoi, cldn       , &
         365             :        state%zm , state%zi    , state%pmid , state%pint , state%u        , state%v , &
         366             :        cam_in%wsx, cam_in%wsy , cam_in%shf , cam_in%cflx(:,1) , wstarent           , &
         367             :        rrho     , ustar       , pblh       , kvm_in     , kvh_in         , kvm     , &
         368             :        kvh      , kvq         , cgh        ,                                         &
         369             :        cgs      , tpert       , qpert      , tke            ,                        &
         370             :        sprod    , sfi         ,                                                      &
         371           0 :        tauresx  , tauresy     , ksrftms    , dragblj )
         372             :
         373             :   ! The diffusivities from diag_TKE can be much larger than from HB in the free
         374             :   ! troposphere and upper atmosphere. These seem to be larger than observations,
         375             :   ! and in WACCM the gw_drag code is already applying an eddy diffusivity in the
         376             :   ! upper atmosphere. Optionally, adjust the diffusivities in the free troposphere
         377             :   ! or the upper atmosphere.
         378             :   !
         379             :   ! NOTE: Further investigation should be done as to why the diffusivities are
         380             :   ! larger in diag_TKE.
         381           0 :   if ((kv_freetrop_scale /= 1._r8) .or. ((kv_top_scale /= 1._r8) .and. (kv_top_pressure > 0._r8))) then
         382           0 :      do i = 1, state%ncol
         383           0 :         do k = 1, pverp
         384             :            ! Outside of the boundary layer?
         385           0 :            if (state%zi(i,k) > pblh(i)) then
         386             :               ! In the upper atmosphere?
         387           0 :               if (state%pint(i,k) <= kv_top_pressure) then
         388           0 :                  kvh(i,k) = kvh(i,k) * kv_top_scale
         389           0 :                  kvm(i,k) = kvm(i,k) * kv_top_scale
         390           0 :                  kvq(i,k) = kvq(i,k) * kv_top_scale
         391             :               else
         392           0 :                  kvh(i,k) = kvh(i,k) * kv_freetrop_scale
         393           0 :                  kvm(i,k) = kvm(i,k) * kv_freetrop_scale
         394           0 :                  kvq(i,k) = kvq(i,k) * kv_freetrop_scale
         395             :               end if
         396             :            else
         397             :               exit
         398             :            end if
         399             :         end do
         400             :      end do
         401             :   end if
         402             :
         403           0 : end subroutine eddy_diff_tend
         404             :
         405             : !=============================================================================== !
         406             : !                                                                                !
         407             : !=============================================================================== !
         408             :
         409           0 : subroutine compute_eddy_diff( pbuf, lchnk  ,                                                      &
         410           0 :                               pcols  , pver   , ncol     , t       , tint, qv       , ztodt   ,   &
         411           0 :                               ql     , qi     , s        , p       , rhoi, cldn     ,             &
         412           0 :                               z      , zi     , pmid     , pi      , u        , v       ,         &
         413           0 :                               taux   , tauy   , shflx    , qflx    , wstarent ,           rrho  , &
         414           0 :                               ustar  , pblh   , kvm_in   , kvh_in  , kvm_out  , kvh_out , kvq   , &
         415           0 :                               cgh    , cgs    , tpert    , qpert   , tke     ,                    &
         416           0 :                               sprod  , sfi    ,                                                   &
         417           0 :                               tauresx, tauresy, ksrftms, dragblj )
         418             :
         419             :   !-------------------------------------------------------------------- !
         420             :   ! Purpose: Interface to compute eddy diffusivities.                   !
         421             :   !          Eddy diffusivities are calculated in a fully implicit way  !
         422             :   !          through iteration process.                                 !
         423             :   ! Author:  Sungsu Park. August. 2006.                                 !
         424             :   !                       May.    2008.                                 !
         425             :   !-------------------------------------------------------------------- !
         426             :
         427           0 :   use diffusion_solver, only: compute_vdiff
         428             :   use cam_history,      only: outfld
         429             :   use phys_debug_util,  only: phys_debug_col
         430             :   use air_composition,  only: cpairv
         431             :   use pbl_utils,        only: calc_ustar, austausch_atm
         432             :   use error_messages,   only: handle_errmsg
         433             :   use coords_1d,        only: Coords1D
         434             :   use wv_saturation,    only: qsat
         435             :   use eddy_diff,        only: trbintd, caleddy
         436             :   use physics_buffer,   only: pbuf_get_field
         437             :
         438             :   ! --------------- !
         439             :   ! Input Variables !
         440             :   ! --------------- !
         441             :
         442             :   type(physics_buffer_desc), pointer, intent(in) :: pbuf(:)
         443             :   integer,  intent(in)    :: lchnk
         444             :   integer,  intent(in)    :: pcols                     ! Number of atmospheric columns [ # ]
         445             :   integer,  intent(in)    :: pver                      ! Number of atmospheric layers  [ # ]
         446             :   integer,  intent(in)    :: ncol                      ! Number of atmospheric columns [ # ]
         447             :   logical,  intent(in)    :: wstarent                  ! .true. means use the 'wstar' entrainment closure.
         448             :   real(r8), intent(in)    :: ztodt                     ! Physics integration time step 2 delta-t [ s ]
         449             :   real(r8), intent(in)    :: t(pcols,pver)             ! Temperature [ K ]
         450             :   real(r8), intent(in)    :: tint(pcols,pver+1)        ! Temperature defined on interfaces [ K ]
         451             :   real(r8), intent(in)    :: qv(pcols,pver)            ! Water vapor  specific humidity [ kg/kg ]
         452             :   real(r8), intent(in)    :: ql(pcols,pver)            ! Liquid water specific humidity [ kg/kg ]
         453             :   real(r8), intent(in)    :: qi(pcols,pver)            ! Ice specific humidity [ kg/kg ]
         454             :   real(r8), intent(in)    :: s(pcols,pver)             ! Dry static energy [ J/kg ]
         455             :   type(Coords1D), intent(in) :: p                      ! Pressure coordinates for solver [ Pa ]
         456             :   real(r8), intent(in)    :: rhoi(pcols,pver+1)        ! Density at interfaces [ kg/m3 ]
         457             :   real(r8), intent(in)    :: cldn(pcols,pver)          ! Stratiform cloud fraction [ fraction ]
         458             :   real(r8), intent(in)    :: z(pcols,pver)             ! Layer mid-point height above surface [ m ]
         459             :   real(r8), intent(in)    :: zi(pcols,pver+1)          ! Interface height above surface [ m ]
         460             :   real(r8), intent(in)    :: pmid(pcols,pver)          ! Layer mid-point pressure [ Pa ]
         461             :   real(r8), intent(in)    :: pi(pcols,pver+1)          ! Interface pressure [ Pa ]
         462             :   real(r8), intent(in)    :: u(pcols,pver)             ! Zonal velocity [ m/s ]
         463             :   real(r8), intent(in)    :: v(pcols,pver)             ! Meridional velocity [ m/s ]
         464             :   real(r8), intent(in)    :: taux(pcols)               ! Zonal wind stress at surface [ N/m2 ]
         465             :   real(r8), intent(in)    :: tauy(pcols)               ! Meridional wind stress at surface [ N/m2 ]
         466             :   real(r8), intent(in)    :: shflx(pcols)              ! Sensible heat flux at surface [ unit ? ]
         467             :   real(r8), intent(in)    :: qflx(pcols)               ! Water vapor flux at surface [ unit ? ]
         468             :   real(r8), intent(in)    :: kvm_in(pcols,pver+1)      ! kvm saved from last timestep [ m2/s ]
         469             :   real(r8), intent(in)    :: kvh_in(pcols,pver+1)      ! kvh saved from last timestep [ m2/s ]
         470             :   real(r8), intent(in)    :: ksrftms(pcols)            ! Surface drag coefficient of turbulent mountain stress [ unit ? ]
         471             :   real(r8), intent(in)    :: dragblj(pcols,pver)       ! Drag profile from Beljaars SGO form drag [ 1/s ]
         472             :
         473             :   ! ---------------- !
         474             :   ! Output Variables !
         475             :   ! ---------------- !
         476             :
         477             :   real(r8), intent(out)   :: kvm_out(pcols,pver+1)     ! Eddy diffusivity for momentum [ m2/s ]
         478             :   real(r8), intent(out)   :: kvh_out(pcols,pver+1)     ! Eddy diffusivity for heat [ m2/s ]
         479             :   real(r8), intent(out)   :: kvq(pcols,pver+1)         ! Eddy diffusivity for constituents, moisture and tracers [ m2/s ]
         480             :                                                        ! (note not having '_out')
         481             :   real(r8), intent(out)   :: rrho(pcols)               ! Reciprocal of density at the lowest layer
         482             :   real(r8), intent(out)   :: ustar(pcols)              ! Surface friction velocity [ m/s ]
         483             :   real(r8), intent(out)   :: pblh(pcols)               ! PBL top height [ m ]
         484             :   real(r8), intent(out)   :: cgh(pcols,pver+1)         ! Counter-gradient term for heat [ J/kg/m ]
         485             :   real(r8), intent(out)   :: cgs(pcols,pver+1)         ! Counter-gradient star [ cg/flux ]
         486             :   real(r8), intent(out)   :: tpert(pcols)              ! Convective temperature excess [ K ]
         487             :   real(r8), intent(out)   :: qpert(pcols)              ! Convective humidity excess [ kg/kg ]
         488             :   real(r8), intent(out)   :: tke(pcols,pver+1)         ! Turbulent kinetic energy [ m2/s2 ]
         489             :   real(r8), intent(out)   :: sprod(pcols,pver+1)       ! Shear production [ m2/s3 ]
         490             :   real(r8), intent(out)   :: sfi(pcols,pver+1)         ! Interfacial layer saturation fraction [ fraction ]
         491             :
         492             :   ! ---------------------- !
         493             :   ! Input-Output Variables !
         494             :   ! ---------------------- !
         495             :
         496             :   real(r8), intent(inout) :: tauresx(pcols)            ! Residual stress to be added in vdiff to correct for turb
         497             :   real(r8), intent(inout) :: tauresy(pcols)            ! Stress mismatch between sfc and atm accumulated in prior timesteps
         498             :
         499             :   ! -------------- !
         500             :   ! pbuf Variables !
         501             :   ! -------------- !
         502             :
         503           0 :   real(r8), pointer :: qrl(:,:)                        ! LW radiative cooling rate
         504           0 :   real(r8), pointer :: wsedl(:,:)                      ! Sedimentation velocity
         505             :                                                        ! of stratiform liquid cloud droplet [ m/s ]
         506             :
         507           0 :   real(r8), pointer :: bprod(:,:)                      ! Buoyancy production of tke [ m2/s3 ]
         508           0 :   integer(i4), pointer :: ipbl(:)                      ! If 1, PBL is CL, while if 0, PBL is STL.
         509           0 :   integer(i4), pointer :: kpblh(:)                     ! Layer index containing PBL top within or at the base interface
         510           0 :   real(r8), pointer :: wstarPBL(:)                     ! Convective velocity within PBL [ m/s ]
         511           0 :   real(r8), pointer :: tkes(:)                         ! TKE at surface interface [ m2/s2 ]
         512           0 :   real(r8), pointer :: went(:)                         ! Entrainment rate at the PBL top interface [ m/s ]
         513             :
         514             :   ! --------------- !
         515             :   ! Local Variables !
         516             :   ! --------------- !
         517             :
         518             :   integer                    icol
         519             :   integer                    i, k, iturb, status
         520             :
         521             :   character(2048)         :: warnstring                ! Warning(s) to print
         522             :   character(128)          :: errstring                 ! Error message
         523             :
         524           0 :   real(r8)                :: kvf(pcols,pver+1)         ! Free atmospheric eddy diffusivity [ m2/s ]
         525           0 :   real(r8)                :: kvm(pcols,pver+1)         ! Eddy diffusivity for momentum [ m2/s ]
         526           0 :   real(r8)                :: kvh(pcols,pver+1)         ! Eddy diffusivity for heat [ m2/s ]
         527             :   real(r8)                :: kvm_preo(pcols,pver+1)    ! Eddy diffusivity for momentum [ m2/s ]
         528             :   real(r8)                :: kvh_preo(pcols,pver+1)    ! Eddy diffusivity for heat [ m2/s ]
         529             :   real(r8)                :: kvm_pre(pcols,pver+1)     ! Eddy diffusivity for momentum [ m2/s ]
         530             :   real(r8)                :: kvh_pre(pcols,pver+1)     ! Eddy diffusivity for heat [ m2/s ]
         531           0 :   real(r8)                :: errorPBL(pcols)           ! Error function showing whether PBL produced convergent solution or not.
         532             :                                                        ! [ unit ? ]
         533           0 :   real(r8)                :: s2(pcols,pver)            ! Shear squared, defined at interfaces except surface [ s-2 ]
         534           0 :   real(r8)                :: n2(pcols,pver)            ! Buoyancy frequency, defined at interfaces except surface [ s-2 ]
         535           0 :   real(r8)                :: ri(pcols,pver)            ! Richardson number, 'n2/s2', defined at interfaces except surface [ s-2 ]
         536           0 :   real(r8)                :: pblhp(pcols)              ! PBL top pressure [ Pa ]
         537           0 :   real(r8)                :: minpblh(pcols)            ! Minimum PBL height based on surface stress
         538             :
         539           0 :   real(r8)                :: qt(pcols,pver)            ! Total specific humidity [ kg/kg ]
         540           0 :   real(r8)                :: sfuh(pcols,pver)          ! Saturation fraction in upper half-layer [ fraction ]
         541           0 :   real(r8)                :: sflh(pcols,pver)          ! Saturation fraction in lower half-layer [ fraction ]
         542           0 :   real(r8)                :: sl(pcols,pver)            ! Liquid water static energy [ J/kg ]
         543           0 :   real(r8)                :: slv(pcols,pver)           ! Liquid water virtual static energy [ J/kg ]
         544           0 :   real(r8)                :: slslope(pcols,pver)       ! Slope of 'sl' in each layer
         545           0 :   real(r8)                :: qtslope(pcols,pver)       ! Slope of 'qt' in each layer
         546           0 :   real(r8)                :: qvfd(pcols,pver)          ! Specific humidity for diffusion [ kg/kg ]
         547           0 :   real(r8)                :: tfd(pcols,pver)           ! Temperature for diffusion [ K ]
         548           0 :   real(r8)                :: slfd(pcols,pver)          ! Liquid static energy [ J/kg ]
         549           0 :   real(r8)                :: qtfd(pcols,pver)          ! Total specific humidity [ kg/kg ]
         550           0 :   real(r8)                :: qlfd(pcols,pver)          ! Liquid water specific humidity for diffusion [ kg/kg ]
         551           0 :   real(r8)                :: ufd(pcols,pver)           ! U-wind for diffusion [ m/s ]
         552           0 :   real(r8)                :: vfd(pcols,pver)           ! V-wind for diffusion [ m/s ]
         553             :
         554             :   ! Buoyancy coefficients : w'b' = ch * w'sl' + cm * w'qt'
         555             :
         556           0 :   real(r8)                :: chu(pcols,pver+1)         ! Heat buoyancy coef for dry states, defined at each interface, finally.
         557           0 :   real(r8)                :: chs(pcols,pver+1)         ! Heat buoyancy coef for sat states, defined at each interface, finally.
         558           0 :   real(r8)                :: cmu(pcols,pver+1)         ! Moisture buoyancy coef for dry states,
         559             :                                                        ! defined at each interface, finally.
         560           0 :   real(r8)                :: cms(pcols,pver+1)         ! Moisture buoyancy coef for sat states,
         561             :                                                        ! defined at each interface, finally.
         562             :
         563           0 :   real(r8)                :: jnk1d(pcols)
         564           0 :   real(r8)                :: jnk2d(pcols,pver+1)
         565           0 :   real(r8)                :: zero(pcols)
         566           0 :   real(r8)                :: zero2d(pcols,pver+1)
         567             :   real(r8)                :: es                     ! Saturation vapor pressure
         568             :   real(r8)                :: qs                     ! Saturation specific humidity
         569             :   real(r8)                :: ep2, templ, temps
         570             :
         571             :   ! ------------------------------- !
         572             :   ! Variables for diagnostic output !
         573             :   ! ------------------------------- !
         574             :
         575           0 :   real(r8)                :: wpert(pcols)              ! Turbulent velocity excess [ m/s ]
         576             :
         577           0 :   real(r8)                :: kbase_o(pcols,ncvmax)     ! Original external base interface index of CL from 'exacol'
         578           0 :   real(r8)                :: ktop_o(pcols,ncvmax)      ! Original external top  interface index of CL from 'exacol'
         579           0 :   real(r8)                :: ncvfin_o(pcols)           ! Original number of CLs from 'exacol'
         580           0 :   real(r8)                :: kbase_mg(pcols,ncvmax)    ! 'kbase' after extending-merging from 'zisocl'
         581           0 :   real(r8)                :: ktop_mg(pcols,ncvmax)     ! 'ktop' after extending-merging from 'zisocl'
         582           0 :   real(r8)                :: ncvfin_mg(pcols)          ! 'ncvfin' after extending-merging from 'zisocl'
         583           0 :   real(r8)                :: kbase_f(pcols,ncvmax)     ! Final 'kbase' after extending-merging & including SRCL
         584           0 :   real(r8)                :: ktop_f(pcols,ncvmax)      ! Final 'ktop' after extending-merging & including SRCL
         585           0 :   real(r8)                :: ncvfin_f(pcols)           ! Final 'ncvfin' after extending-merging & including SRCL
         586           0 :   real(r8)                :: wet(pcols,ncvmax)         ! Entrainment rate at the CL top  [ m/s ]
         587           0 :   real(r8)                :: web(pcols,ncvmax)         ! Entrainment rate at the CL base [ m/s ].
         588             :                                                        ! Set to zero if CL is based at surface.
         589           0 :   real(r8)                :: jtbu(pcols,ncvmax)        ! Buoyancy jump across the CL top  [ m/s2 ]
         590           0 :   real(r8)                :: jbbu(pcols,ncvmax)        ! Buoyancy jump across the CL base [ m/s2 ]
         591           0 :   real(r8)                :: evhc(pcols,ncvmax)        ! Evaporative enhancement factor at the CL top
         592           0 :   real(r8)                :: jt2slv(pcols,ncvmax)      ! Jump of slv ( across two layers ) at CL top used only for evhc [ J/kg ]
         593           0 :   real(r8)                :: n2ht(pcols,ncvmax)        ! n2 defined at the CL top  interface but using
         594             :                                                        ! sfuh(kt)   instead of sfi(kt) [ s-2 ]
         595           0 :   real(r8)                :: n2hb(pcols,ncvmax)        ! n2 defined at the CL base interface but using
         596             :                                                        ! sflh(kb-1) instead of sfi(kb) [ s-2 ]
         597           0 :   real(r8)                :: lwp(pcols,ncvmax)         ! LWP in the CL top layer [ kg/m2 ]
         598           0 :   real(r8)                :: opt_depth(pcols,ncvmax)   ! Optical depth of the CL top layer
         599           0 :   real(r8)                :: radinvfrac(pcols,ncvmax)  ! Fraction of radiative cooling confined in the top portion of CL top layer
         600           0 :   real(r8)                :: radf(pcols,ncvmax)        ! Buoyancy production at the CL top due to LW radiative cooling [ m2/s3 ]
         601           0 :   real(r8)                :: wstar(pcols,ncvmax)       ! Convective velocity in each CL [ m/s ]
         602           0 :   real(r8)                :: wstar3fact(pcols,ncvmax)  ! Enhancement of 'wstar3' due to entrainment (inverse) [ no unit ]
         603           0 :   real(r8)                :: ebrk(pcols,ncvmax)        ! Net mean TKE of CL including entrainment effect [ m2/s2 ]
         604           0 :   real(r8)                :: wbrk(pcols,ncvmax)        ! Net mean normalized TKE (W) of CL,
         605             :                                                        ! 'ebrk/b1' including entrainment effect [ m2/s2 ]
         606           0 :   real(r8)                :: lbrk(pcols,ncvmax)        ! Energetic internal thickness of CL [m]
         607           0 :   real(r8)                :: ricl(pcols,ncvmax)        ! CL internal mean Richardson number
         608           0 :   real(r8)                :: ghcl(pcols,ncvmax)        ! Half of normalized buoyancy production of CL
         609           0 :   real(r8)                :: shcl(pcols,ncvmax)        ! Galperin instability function of heat-moisture of CL
         610           0 :   real(r8)                :: smcl(pcols,ncvmax)        ! Galperin instability function of mementum of CL
         611           0 :   real(r8)                :: ghi(pcols,pver+1)         ! Half of normalized buoyancy production at all interfaces
         612           0 :   real(r8)                :: shi(pcols,pver+1)         ! Galperin instability function of heat-moisture at all interfaces
         613           0 :   real(r8)                :: smi(pcols,pver+1)         ! Galperin instability function of heat-moisture at all interfaces
         614           0 :   real(r8)                :: rii(pcols,pver+1)         ! Interfacial Richardson number defined at all interfaces
         615           0 :   real(r8)                :: lengi(pcols,pver+1)       ! Turbulence length scale at all interfaces [ m ]
         616           0 :   real(r8)                :: wcap(pcols,pver+1)        ! Normalized TKE at all interfaces [ m2/s2 ]
         617             :   ! For sedimentation-entrainment feedback
         618           0 :   real(r8)                :: wsed(pcols,ncvmax)        ! Sedimentation velocity at the top of each CL [ m/s ]
         619             :
         620           0 :   integer(i4)             :: turbtype(pcols,pver+1)    ! Turbulence type identifier at all interfaces [ no unit ]
         621             :
         622             :   ! ---------- !
         623             :   ! Parameters !
         624             :   ! ---------- !
         625             :
         626             :   logical,          parameter :: use_kvf        =  .false.      ! .true. (.false.) : initialize kvh/kvm =  kvf ( 0. )
         627             :   real(r8),         parameter :: lambda         =   0.5_r8      ! Under-relaxation factor ( 0 < lambda =< 1 )
         628             :
         629             :   ! ---------- !
         630             :   ! Initialize !
         631             :   ! ---------- !
         632             :
         633           0 :   zero(:)     = 0._r8
         634           0 :   zero2d(:,:) = 0._r8
         635             :
         636             :   ! ---------------------------------------------- !
         637             :   ! Get LW radiative heating out of physics buffer !
         638             :   ! ---------------------------------------------- !
         639           0 :   call pbuf_get_field(pbuf, qrl_idx,   qrl)
         640           0 :   call pbuf_get_field(pbuf, wsedl_idx, wsedl)
         641             :
         642             :   ! These fields are put into the pbuf for UNICON only.
         643           0 :   if (unicon_is_on) then
         644           0 :      call pbuf_get_field(pbuf, bprod_idx,    bprod)
         645           0 :      call pbuf_get_field(pbuf, ipbl_idx,     ipbl)
         646           0 :      call pbuf_get_field(pbuf, kpblh_idx,    kpblh)
         647           0 :      call pbuf_get_field(pbuf, wstarPBL_idx, wstarPBL)
         648           0 :      call pbuf_get_field(pbuf, tkes_idx,     tkes)
         649           0 :      call pbuf_get_field(pbuf, went_idx,     went)
         650             :   else
         651           0 :      allocate(bprod(pcols,pverp), ipbl(pcols), kpblh(pcols), wstarPBL(pcols), tkes(pcols), went(pcols))
         652             :   end if
         653             :
         654             :   ! ----------------------- !
         655             :   ! Main Computation Begins !
         656             :   ! ----------------------- !
         657             :
         658           0 :   ufd(:ncol,:)  = u(:ncol,:)
         659           0 :   vfd(:ncol,:)  = v(:ncol,:)
         660           0 :   tfd(:ncol,:)  = t(:ncol,:)
         661           0 :   qvfd(:ncol,:) = qv(:ncol,:)
         662           0 :   qlfd(:ncol,:) = ql(:ncol,:)
         663             :
         664           0 :   do iturb = 1, nturb
         665             :
         666             :      ! Total stress includes 'tms'.
         667             :      ! Here, in computing 'tms', we can use either iteratively changed 'ufd,vfd' or the
         668             :      ! initially given 'u,v' to the PBL scheme. Note that normal stress, 'taux, tauy'
         669             :      ! are not changed by iteration. In order to treat 'tms' in a fully implicit way,
         670             :      ! I am using updated wind, here.
         671             :
         672             :      ! Compute ustar
         673             :      call calc_ustar( ncol, tfd(:ncol,pver), pmid(:ncol,pver), &
         674           0 :                       taux(:ncol) - ksrftms(:ncol) * ufd(:ncol,pver), & ! Zonal wind stress
         675             :                       tauy(:ncol) - ksrftms(:ncol) * vfd(:ncol,pver), & ! Meridional wind stress
         676           0 :                       rrho(:ncol), ustar(:ncol))
         677           0 :      minpblh(:ncol) = 100.0_r8 * ustar(:ncol)   ! By construction, 'minpblh' is larger than 1 [m] when 'ustar_min = 0.01'.
         678             :
         679             :      ! Calculate (qt,sl,n2,s2,ri) from a given set of (t,qv,ql,qi,u,v)
         680             :
         681             :      call trbintd( &
         682             :                    pcols    , pver    , ncol  , z       , ufd     , vfd     , tfd   , pmid    , &
         683             :                    s2       , n2      , ri    , zi      , pi      , cldn    , qtfd  , qvfd    , &
         684             :                    qlfd     , qi      , sfi   , sfuh    , sflh    , slfd    , slv   , slslope , &
         685           0 :                    qtslope  , chs     , chu   , cms     , cmu     )
         686             :
         687             :      ! Save initial (i.e., before iterative diffusion) profile of (qt,sl) at each iteration.
         688             :      ! Only necessary for (qt,sl) not (u,v) because (qt,sl) are newly calculated variables.
         689             :
         690           0 :      if( iturb == 1 ) then
         691           0 :         qt(:ncol,:) = qtfd(:ncol,:)
         692           0 :         sl(:ncol,:) = slfd(:ncol,:)
         693             :      endif
         694             :
         695             :      ! Get free atmosphere exchange coefficients. This 'kvf' is not used in UW moist PBL scheme
         696             :      if (use_kvf) then
         697             :         call austausch_atm(pcols, ncol, pver, ntop_eddy, nbot_eddy, &
         698             :              ml2, ri, s2, kvf )
         699             :      else
         700           0 :         kvf = 0._r8
         701             :      end if
         702             :
         703             :      ! Initialize kvh/kvm to send to caleddy, depending on model timestep and iteration number
         704             :      ! This is necessary for 'wstar-based' entrainment closure.
         705             :
         706           0 :      if( iturb == 1 ) then
         707           0 :         if( is_first_step() ) then
         708             :            ! First iteration of first model timestep : Use free tropospheric value or zero.
         709           0 :            kvh(:ncol,:) = kvf(:ncol,:)
         710           0 :            kvm(:ncol,:) = kvf(:ncol,:)
         711             :         else
         712             :            ! First iteration on any model timestep except the first : Use value from previous timestep
         713           0 :            kvh(:ncol,:) = kvh_in(:ncol,:)
         714           0 :            kvm(:ncol,:) = kvm_in(:ncol,:)
         715             :         endif
         716             :      else
         717             :         ! Not the first iteration : Use from previous iteration
         718           0 :         kvh(:ncol,:) = kvh_out(:ncol,:)
         719           0 :         kvm(:ncol,:) = kvm_out(:ncol,:)
         720             :      endif
         721             :
         722             :      ! Calculate eddy diffusivity (kvh_out,kvm_out) and (tke,bprod,sprod) using
         723             :      ! a given (kvh,kvm) which are used only for initializing (bprod,sprod)  at
         724             :      ! the first part of caleddy. (bprod,sprod) are fully updated at the end of
         725             :      ! caleddy after calculating (kvh_out,kvm_out)
         726             :
         727             :      call caleddy( pcols     , pver      , ncol      ,                     &
         728             :                    slfd      , qtfd      , qlfd      , slv      ,ufd     , &
         729             :                    vfd       , pi        , z         , zi       ,          &
         730             :                    qflx      , shflx     , slslope   , qtslope  ,          &
         731             :                    chu       , chs       , cmu       , cms      ,sfuh    , &
         732             :                    sflh      , n2        , s2        , ri       ,rrho    , &
         733             :                    pblh      , ustar     ,                                 &
         734             :                    kvh       , kvm       , kvh_out   , kvm_out  ,          &
         735             :                    tpert     , qpert     , qrl       , kvf      , tke    , &
         736             :                    wstarent  , bprod     , sprod     , minpblh  , wpert  , &
         737             :                    tkes      , went      , turbtype  ,                     &
         738             :                    kbase_o   , ktop_o    , ncvfin_o  ,                     &
         739             :                    kbase_mg  , ktop_mg   , ncvfin_mg ,                     &
         740             :                    kbase_f   , ktop_f    , ncvfin_f  ,                     &
         741             :                    wet       , web       , jtbu      , jbbu     ,          &
         742             :                    evhc      , jt2slv    , n2ht      , n2hb     ,          &
         743             :                    lwp       , opt_depth , radinvfrac, radf     ,          &
         744             :                    wstar     , wstar3fact,                                 &
         745             :                    ebrk      , wbrk      , lbrk      , ricl     , ghcl   , &
         746             :                    shcl      , smcl      , ghi       , shi      , smi    , &
         747             :                    rii       , lengi     , wcap      , pblhp    , cldn   , &
         748             :                    ipbl      , kpblh     , wsedl     , wsed, &
         749           0 :                    warnstring, errstring)
         750             :
         751           0 :      if (trim(warnstring) /= "") then
         752           0 :         write(iulog,*) "eddy_diff_cam: Messages from caleddy follow."
         753           0 :         write(iulog,*) warnstring
         754             :      end if
         755             :
         756           0 :      call handle_errmsg(errstring, subname="caleddy")
         757             :
         758             :      ! Calculate errorPBL to check whether PBL produced convergent solutions or not.
         759             :
         760           0 :      if( iturb == nturb ) then
         761           0 :         do i = 1, ncol
         762           0 :            errorPBL(i) = 0._r8
         763           0 :            do k = 1, pver
         764           0 :               errorPBL(i) = errorPBL(i) + ( kvh(i,k) - kvh_out(i,k) )**2
         765             :            end do
         766           0 :            errorPBL(i) = sqrt(errorPBL(i)/pver)
         767             :         end do
         768             :      end if
         769             :
         770             :      ! Eddy diffusivities which will be used for the initialization of (bprod,
         771             :      ! sprod) in 'caleddy' at the next iteration step.
         772             :
         773           0 :      if( iturb > 1 .and. iturb < nturb ) then
         774           0 :         kvm_out(:ncol,:) = lambda * kvm_out(:ncol,:) + ( 1._r8 - lambda ) * kvm(:ncol,:)
         775           0 :         kvh_out(:ncol,:) = lambda * kvh_out(:ncol,:) + ( 1._r8 - lambda ) * kvh(:ncol,:)
         776             :      endif
         777             :
         778             :      ! Set nonlocal terms to zero for flux diagnostics, since not used by caleddy.
         779             :
         780           0 :      cgh(:ncol,:) = 0._r8
         781           0 :      cgs(:ncol,:) = 0._r8
         782             :
         783           0 :      if( iturb < nturb ) then
         784             :
         785             :         ! Each time we diffuse the original state
         786             :
         787           0 :         slfd(:ncol,:)  = sl(:ncol,:)
         788           0 :         qtfd(:ncol,:)  = qt(:ncol,:)
         789           0 :         ufd(:ncol,:)   = u(:ncol,:)
         790           0 :         vfd(:ncol,:)   = v(:ncol,:)
         791             :
         792             :         ! Diffuse initial profile of each time step using a given (kvh_out,kvm_out)
         793             :         ! In the below 'compute_vdiff', (slfd,qtfd,ufd,vfd) are 'inout' variables.
         794             :
         795             :         call compute_vdiff( lchnk   ,                                                  &
         796             :                             pcols   , pver     , 1        , ncol         , tint, &
         797             :                             p        , t        , rhoi, ztodt        , taux      , &
         798             :                             tauy    , shflx    , qflx     , &
         799             :                             kvh_out , kvm_out  , kvh_out  , cgs          , cgh       , &
         800             :                             zi      , ksrftms  , dragblj  , &
         801             :                             zero    , fieldlist_wet, fieldlist_molec, &
         802             :                             ufd     , vfd      , qtfd     , slfd         ,             &
         803             :                             jnk1d   , jnk1d    , jnk2d    , jnk1d        , errstring , &
         804             :                             tauresx , tauresy  , 0        , cpairv(:,:,lchnk), zero, &
         805           0 :                             .false., .false. )
         806             :
         807             :         call handle_errmsg(errstring, subname="compute_vdiff", &
         808           0 :              extra_msg="compute_vdiff called from eddy_diff_cam")
         809             :
         810             :         ! Retrieve (tfd,qvfd,qlfd) from (slfd,qtfd) in order to
         811             :         ! use 'trbintd' at the next iteration.
         812             :
         813           0 :         do k = 1, pver
         814           0 :            do i = 1, ncol
         815             :               ! ----------------------------------------------------- !
         816             :               ! Compute the condensate 'qlfd' in the updated profiles !
         817             :               ! ----------------------------------------------------- !
         818             :               ! Option.1 : Assume grid-mean condensate is homogeneously diffused by the moist turbulence scheme.
         819             :               !            This should be used if 'pseudodiff = .false.' in vertical_diffusion.F90.
         820             :               ! Modification : Need to be check whether below is correct in the presence of ice, qi.
         821             :               !                I should understand why the variation of ice, qi is neglected during diffusion.
         822           0 :               templ     = ( slfd(i,k) - gravit*z(i,k) ) / cpair
         823           0 :               call qsat( templ, pmid(i,k), es, qs)
         824           0 :               ep2       =  .622_r8
         825           0 :               temps     =   templ + ( qtfd(i,k) - qs ) / ( cpair / latvap + latvap * qs / ( rair * templ**2 ) )
         826           0 :               call qsat( temps, pmid(i,k), es, qs)
         827           0 :               qlfd(i,k) =   max( qtfd(i,k) - qi(i,k) - qs ,0._r8 )
         828             :               ! Option.2 : Assume condensate is not diffused by the moist turbulence scheme.
         829             :               !            This should bs used if 'pseudodiff = .true.'  in vertical_diffusion.F90.
         830             :               ! qlfd(i,k) = ql(i,k)
         831             :               ! ----------------------------- !
         832             :               ! Compute the other 'qvfd, tfd' !
         833             :               ! ----------------------------- !
         834           0 :               qvfd(i,k) = max( 0._r8, qtfd(i,k) - qi(i,k) - qlfd(i,k) )
         835           0 :               tfd(i,k)  = ( slfd(i,k) + latvap * qlfd(i,k) + (latvap+latice) * qi(i,k) - gravit*z(i,k)) / cpair
         836             :            end do
         837             :         end do
         838             :      endif
         839             :
         840             :   end do  ! End of 'iturb' iteration
         841             :
         842           0 :   kvq(:ncol,:) = kvh_out(:ncol,:)
         843             :
         844             :   ! Compute 'wstar' within the PBL for use in the future convection scheme.
         845             :
         846           0 :   do i = 1, ncol
         847           0 :      if(ipbl(i) == 1) then
         848           0 :         wstarPBL(i) = max( 0._r8, wstar(i,1) )
         849             :      else
         850           0 :         wstarPBL(i) = 0._r8
         851             :      endif
         852             :   end do
         853             :
         854             :   ! --------------------------------------------------------------- !
         855             :   ! Writing for detailed diagnostic analysis of UW moist PBL scheme !
         856             :   ! --------------------------------------------------------------- !
         857             :
         858           0 :   call outfld( 'WGUSTD' , wpert, pcols, lchnk )
         859             :
         860           0 :   call outfld( 'BPROD   ', bprod, pcols, lchnk )
         861           0 :   call outfld( 'SPROD   ', sprod, pcols, lchnk )
         862           0 :   call outfld( 'SFI     ', sfi,   pcols, lchnk )
         863             :
         864           0 :   call outfld( 'UW_errorPBL',    errorPBL,   pcols,   lchnk )
         865             :
         866           0 :   call outfld( 'UW_n2',          n2,         pcols,   lchnk )
         867           0 :   call outfld( 'UW_s2',          s2,         pcols,   lchnk )
         868           0 :   call outfld( 'UW_ri',          ri,         pcols,   lchnk )
         869             :
         870           0 :   call outfld( 'UW_sfuh',        sfuh,       pcols,   lchnk )
         871           0 :   call outfld( 'UW_sflh',        sflh,       pcols,   lchnk )
         872           0 :   call outfld( 'UW_sfi',         sfi,        pcols,   lchnk )
         873             :
         874           0 :   call outfld( 'UW_cldn',        cldn,       pcols,   lchnk )
         875           0 :   call outfld( 'UW_qrl',         qrl,        pcols,   lchnk )
         876           0 :   call outfld( 'UW_ql',          qlfd,       pcols,   lchnk )
         877             :
         878           0 :   call outfld( 'UW_chu',         chu,        pcols,   lchnk )
         879           0 :   call outfld( 'UW_chs',         chs,        pcols,   lchnk )
         880           0 :   call outfld( 'UW_cmu',         cmu,        pcols,   lchnk )
         881           0 :   call outfld( 'UW_cms',         cms,        pcols,   lchnk )
         882             :
         883           0 :   call outfld( 'UW_tke',         tke,        pcols,   lchnk )
         884           0 :   call outfld( 'UW_wcap',        wcap,       pcols,   lchnk )
         885           0 :   call outfld( 'UW_bprod',       bprod,      pcols,   lchnk )
         886           0 :   call outfld( 'UW_sprod',       sprod,      pcols,   lchnk )
         887             :
         888           0 :   call outfld( 'UW_kvh',         kvh_out,    pcols,   lchnk )
         889           0 :   call outfld( 'UW_kvm',         kvm_out,    pcols,   lchnk )
         890             :
         891           0 :   call outfld( 'UW_pblh',        pblh,       pcols,   lchnk )
         892           0 :   call outfld( 'UW_pblhp',       pblhp,      pcols,   lchnk )
         893           0 :   call outfld( 'UW_tpert',       tpert,      pcols,   lchnk )
         894           0 :   call outfld( 'UW_qpert',       qpert,      pcols,   lchnk )
         895           0 :   call outfld( 'UW_wpert',       wpert,      pcols,   lchnk )
         896             :
         897           0 :   call outfld( 'UW_ustar',       ustar,      pcols,   lchnk )
         898           0 :   call outfld( 'UW_tkes',        tkes,       pcols,   lchnk )
         899           0 :   call outfld( 'UW_minpblh',     minpblh,    pcols,   lchnk )
         900             :
         901           0 :   call outfld( 'UW_turbtype',    real(turbtype,r8),   pcols,   lchnk )
         902             :
         903           0 :   call outfld( 'UW_kbase_o',     kbase_o,    pcols,   lchnk )
         904           0 :   call outfld( 'UW_ktop_o',      ktop_o,     pcols,   lchnk )
         905           0 :   call outfld( 'UW_ncvfin_o',    ncvfin_o,   pcols,   lchnk )
         906             :
         907           0 :   call outfld( 'UW_kbase_mg',    kbase_mg,   pcols,   lchnk )
         908           0 :   call outfld( 'UW_ktop_mg',     ktop_mg,    pcols,   lchnk )
         909           0 :   call outfld( 'UW_ncvfin_mg',   ncvfin_mg,  pcols,   lchnk )
         910             :
         911           0 :   call outfld( 'UW_kbase_f',     kbase_f,    pcols,   lchnk )
         912           0 :   call outfld( 'UW_ktop_f',      ktop_f,     pcols,   lchnk )
         913           0 :   call outfld( 'UW_ncvfin_f',    ncvfin_f,   pcols,   lchnk )
         914             :
         915           0 :   call outfld( 'UW_wet',         wet,        pcols,   lchnk )
         916           0 :   call outfld( 'UW_web',         web,        pcols,   lchnk )
         917           0 :   call outfld( 'UW_jtbu',        jtbu,       pcols,   lchnk )
         918           0 :   call outfld( 'UW_jbbu',        jbbu,       pcols,   lchnk )
         919           0 :   call outfld( 'UW_evhc',        evhc,       pcols,   lchnk )
         920           0 :   call outfld( 'UW_jt2slv',      jt2slv,     pcols,   lchnk )
         921           0 :   call outfld( 'UW_n2ht',        n2ht,       pcols,   lchnk )
         922           0 :   call outfld( 'UW_n2hb',        n2hb,       pcols,   lchnk )
         923           0 :   call outfld( 'UW_lwp',         lwp,        pcols,   lchnk )
         924           0 :   call outfld( 'UW_optdepth',    opt_depth,  pcols,   lchnk )
         925           0 :   call outfld( 'UW_radfrac',     radinvfrac, pcols,   lchnk )
         926           0 :   call outfld( 'UW_radf',        radf,       pcols,   lchnk )
         927           0 :   call outfld( 'UW_wstar',       wstar,      pcols,   lchnk )
         928           0 :   call outfld( 'UW_wstar3fact',  wstar3fact, pcols,   lchnk )
         929           0 :   call outfld( 'UW_ebrk',        ebrk,       pcols,   lchnk )
         930           0 :   call outfld( 'UW_wbrk',        wbrk,       pcols,   lchnk )
         931           0 :   call outfld( 'UW_lbrk',        lbrk,       pcols,   lchnk )
         932           0 :   call outfld( 'UW_ricl',        ricl,       pcols,   lchnk )
         933           0 :   call outfld( 'UW_ghcl',        ghcl,       pcols,   lchnk )
         934           0 :   call outfld( 'UW_shcl',        shcl,       pcols,   lchnk )
         935           0 :   call outfld( 'UW_smcl',        smcl,       pcols,   lchnk )
         936             :
         937           0 :   call outfld( 'UW_gh',          ghi,        pcols,   lchnk )
         938           0 :   call outfld( 'UW_sh',          shi,        pcols,   lchnk )
         939           0 :   call outfld( 'UW_sm',          smi,        pcols,   lchnk )
         940           0 :   call outfld( 'UW_ria',         rii,        pcols,   lchnk )
         941           0 :   call outfld( 'UW_leng',        lengi,      pcols,   lchnk )
         942             :
         943           0 :   call outfld( 'UW_wsed',        wsed,       pcols,   lchnk )
         944             :
         945           0 :   if (.not. unicon_is_on) then
         946           0 :      deallocate(bprod, ipbl, kpblh, wstarPBL, tkes, went)
         947             :   end if
         948             :
         949           0 : end subroutine compute_eddy_diff
         950             :
         951             : end module eddy_diff_cam

Generated by: LCOV version 1.14