LCOV - code coverage report
Current view: | top level - physics/cam - vertical_diffusion.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 270 | 615 | 43.9 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 5 | 6 | 83.3 %

              Line data    Source code
           1             : module vertical_diffusion
           2             :
           3             : !----------------------------------------------------------------------------------------------------- !
           4             : ! Module to compute vertical diffusion of momentum,  moisture, trace constituents                      !
           5             : ! and static energy. Separate modules compute                                                          !
           6             : !   1. stresses associated with turbulent flow over orography                                          !
           7             : !      ( turbulent mountain stress )                                                                   !
           8             : !   2. eddy diffusivities, including nonlocal tranport terms                                           !
           9             : !   3. molecular diffusivities                                                                         !
          10             : ! Lastly, a implicit diffusion solver is called, and tendencies retrieved by                           !
          11             : ! differencing the diffused and initial states.                                                        !
          12             : !                                                                                                      !
          13             : ! Calling sequence:                                                                                    !
          14             : !                                                                                                      !
          15             : !  vertical_diffusion_init      Initializes vertical diffustion constants and modules                  !
          16             : !        init_molec_diff        Initializes molecular diffusivity module                               !
          17             : !        init_eddy_diff         Initializes eddy diffusivity module (includes PBL)                     !
          18             : !        init_tms               Initializes turbulent mountain stress module                           !
          19             : !        init_vdiff             Initializes diffusion solver module                                    !
          20             : !  vertical_diffusion_ts_init   Time step initialization (only used for upper boundary condition)      !
          21             : !  vertical_diffusion_tend      Computes vertical diffusion tendencies                                 !
          22             : !        compute_tms            Computes turbulent mountain stresses                                   !
          23             : !        compute_eddy_diff      Computes eddy diffusivities and countergradient terms                  !
          24             : !        compute_vdiff          Solves vertical diffusion equations, including molecular diffusivities !
          25             : !                                                                                                      !
          26             : !----------------------------------------------------------------------------------------------------- !
          27             : ! Some notes on refactoring changes made in 2015, which were not quite finished.                       !
          28             : !                                                                                                      !
          29             : !      - eddy_diff_tend should really only have state, pbuf, and cam_in as inputs. The process of      !
          30             : !        removing these arguments, and referring to pbuf fields instead, is not complete.              !
          31             : !                                                                                                      !
          32             : !      - compute_vdiff was intended to be split up into three components:                              !
          33             : !                                                                                                      !
          34             : !         1. Diffusion of winds and heat ("U", "V", and "S" in the fieldlist object).                  !
          35             : !                                                                                                      !
          36             : !         2. Turbulent diffusion of a single constituent                                               !
          37             : !                                                                                                      !
          38             : !         3. Molecular diffusion of a single constituent                                               !
          39             : !                                                                                                      !
          40             : !        This reorganization would allow the three resulting functions to each use a simpler interface !
          41             : !        than the current combined version, and possibly also remove the need to use the fieldlist     !
          42             : !        object at all.                                                                                !
          43             : !                                                                                                      !
          44             : !      - The conditionals controlled by "do_pbl_diags" are somewhat scattered. It might be better to   !
          45             : !        pull out these diagnostic calculations and outfld calls into separate functions.              !
          46             : !                                                                                                      !
          47             : !---------------------------Code history-------------------------------------------------------------- !
          48             : ! J. Rosinski : Jun. 1992                                                                              !
          49             : ! J. McCaa    : Sep. 2004                                                                              !
          50             : ! S. Park     : Aug. 2006, Dec. 2008. Jan. 2010                                                        !
          51             : !----------------------------------------------------------------------------------------------------- !
          52             :
          53             : use shr_kind_mod,     only : r8 => shr_kind_r8, i4=> shr_kind_i4
          54             : use ppgrid,           only : pcols, pver, pverp
          55             : use constituents,     only : pcnst
          56             : use diffusion_solver, only : vdiff_selector
          57             : use cam_abortutils,   only : endrun
          58             : use error_messages,   only : handle_errmsg
          59             : use physconst,        only :          &
          60             :      cpair  , &     ! Specific heat of dry air
          61             :      gravit , &     ! Acceleration due to gravity
          62             :      rair   , &     ! Gas constant for dry air
          63             :      zvir   , &     ! rh2o/rair - 1
          64             :      latvap , &     ! Latent heat of vaporization
          65             :      latice , &     ! Latent heat of fusion
          66             :      karman , &     ! von Karman constant
          67             :      mwdry  , &     ! Molecular weight of dry air
          68             :      avogad         ! Avogadro's number
          69             : use cam_history,      only : fieldname_len
          70             : use perf_mod
          71             : use cam_logfile,      only : iulog
          72             : use ref_pres,         only : do_molec_diff, nbot_molec
          73             : use phys_control,     only : phys_getopts
          74             : use time_manager,     only : is_first_step
          75             : implicit none
          76             : private
          77             : save
          78             :
          79             : ! ----------------- !
          80             : ! Public interfaces !
          81             : ! ----------------- !
          82             :
          83             : public vd_readnl
          84             : public vd_register                                   ! Register multi-time-level variables with physics buffer
          85             : public vertical_diffusion_init                       ! Initialization
          86             : public vertical_diffusion_ts_init                    ! Time step initialization (only used for upper boundary condition)
          87             : public vertical_diffusion_tend                       ! Full vertical diffusion routine
          88             :
          89             : ! ------------ !
          90             : ! Private data !
          91             : ! ------------ !
          92             :
          93             : character(len=16)    :: eddy_scheme                  ! Default set in phys_control.F90, use namelist to change
          94             : !     'HB'       = Holtslag and Boville (default)
          95             : !     'HBR'      = Holtslag and Boville and Rash
          96             : !     'diag_TKE' = Bretherton and Park ( UW Moist Turbulence Scheme )
          97             : logical, parameter   :: wstarent = .true.            ! Use wstar (.true.) or TKE (.false.) entrainment closure
          98             : ! ( when 'diag_TKE' scheme is selected )
          99             : logical              :: do_pseudocon_diff = .false.  ! If .true., do pseudo-conservative variables diffusion
         100             :
         101             : character(len=16)    :: shallow_scheme               ! Shallow convection scheme
         102             :
         103             : type(vdiff_selector) :: fieldlist_wet                ! Logical switches for moist mixing ratio diffusion
         104             : type(vdiff_selector) :: fieldlist_dry                ! Logical switches for dry mixing ratio diffusion
         105             : type(vdiff_selector) :: fieldlist_molec              ! Logical switches for molecular diffusion
         106             : integer              :: tke_idx, kvh_idx, kvm_idx    ! TKE and eddy diffusivity indices for fields in the physics buffer
         107             : integer              :: kvt_idx                      ! Index for kinematic molecular conductivity
         108             : integer              :: tauresx_idx, tauresy_idx     ! Redisual stress for implicit surface stress
         109             :
         110             : character(len=fieldname_len) :: vdiffnam(pcnst)      ! Names of vertical diffusion tendencies
         111             : integer              :: ixcldice, ixcldliq           ! Constituent indices for cloud liquid and ice water
         112             : integer              :: ixnumice, ixnumliq
         113             :
         114             : integer              :: pblh_idx, tpert_idx, qpert_idx
         115             :
         116             : ! pbuf fields for unicon
         117             : integer              :: qtl_flx_idx  = -1            ! for use in cloud macrophysics when UNICON is on
         118             : integer              :: qti_flx_idx  = -1            ! for use in cloud macrophysics when UNICON is on
         119             :
         120             : ! pbuf fields for tms
         121             : integer              :: ksrftms_idx  = -1
         122             : integer              :: tautmsx_idx  = -1
         123             : integer              :: tautmsy_idx  = -1
         124             :
         125             : ! pbuf fields for blj (Beljaars)
         126             : integer              :: dragblj_idx  = -1
         127             : integer              :: taubljx_idx  = -1
         128             : integer              :: taubljy_idx  = -1
         129             :
         130             : ! pbuf field for clubb top above which HB (Holtslag Boville) scheme may be enabled
         131             : integer              :: clubbtop_idx = -1
         132             :
         133             : logical              :: diff_cnsrv_mass_check        ! do mass conservation check
         134             : logical              :: do_iss                       ! switch for implicit turbulent surface stress
         135             : logical              :: prog_modal_aero = .false.    ! set true if prognostic modal aerosols are present
         136             : integer              :: pmam_ncnst = 0               ! number of prognostic modal aerosol constituents
         137             : integer, allocatable :: pmam_cnst_idx(:)             ! constituent indices of prognostic modal aerosols
         138             :
         139             : logical              :: do_pbl_diags = .false.
         140             : logical              :: waccmx_mode = .false.
         141             : logical              :: do_hb_above_clubb = .false.
         142             :
         143             : real(r8),allocatable :: kvm_sponge(:)
         144             :
         145             : contains
         146             :
         147             :   ! =============================================================================== !
         148             :   !                                                                                 !
         149             :   ! =============================================================================== !
         150        1536 : subroutine vd_readnl(nlfile)
         151             :
         152             :   use namelist_utils,  only: find_group_name
         153             :   use units,           only: getunit, freeunit
         154             :   use spmd_utils,      only: masterproc, masterprocid, mpi_logical, mpicom
         155             :   use shr_log_mod,     only: errMsg => shr_log_errMsg
         156             :   use trb_mtn_stress_cam, only: trb_mtn_stress_readnl
         157             :   use beljaars_drag_cam, only: beljaars_drag_readnl
         158             :   use eddy_diff_cam,   only: eddy_diff_readnl
         159             :
         160             :   character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         161             :
         162             :   ! Local variables
         163             :   integer :: unitn, ierr
         164             :   character(len=*), parameter :: subname = 'vd_readnl'
         165             :
         166             :   namelist /vert_diff_nl/ diff_cnsrv_mass_check, do_iss
         167             :   !-----------------------------------------------------------------------------
         168             :
         169        1536 :   if (masterproc) then
         170           2 :      unitn = getunit()
         171           2 :      open( unitn, file=trim(nlfile), status='old' )
         172           2 :      call find_group_name(unitn, 'vert_diff_nl', status=ierr)
         173           2 :      if (ierr == 0) then
         174           2 :         read(unitn, vert_diff_nl, iostat=ierr)
         175           2 :         if (ierr /= 0) then
         176           0 :            call endrun(subname // ':: ERROR reading namelist')
         177             :         end if
         178             :      end if
         179           2 :      close(unitn)
         180           2 :      call freeunit(unitn)
         181             :   end if
         182             :
         183        1536 :   call mpi_bcast(diff_cnsrv_mass_check, 1, mpi_logical, masterprocid, mpicom, ierr)
         184        1536 :   if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
         185        1536 :   call mpi_bcast(do_iss,                1, mpi_logical, masterprocid, mpicom, ierr)
         186        1536 :   if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
         187             :
         188             :   ! Get eddy_scheme setting from phys_control.
         189             :   call phys_getopts( eddy_scheme_out          =          eddy_scheme, &
         190        1536 :        shallow_scheme_out       =       shallow_scheme )
         191             :
         192             :   ! TMS reads its own namelist.
         193        1536 :   call trb_mtn_stress_readnl(nlfile)
         194             :
         195             :   ! Beljaars reads its own namelist.
         196        1536 :   call beljaars_drag_readnl(nlfile)
         197             :
         198        1536 :   if (eddy_scheme == 'diag_TKE' .or. eddy_scheme == 'SPCAM_m2005' ) call eddy_diff_readnl(nlfile)
         199             :
         200        1536 : end subroutine vd_readnl
         201             :
         202             : ! =============================================================================== !
         203             : !                                                                                 !
         204             : ! =============================================================================== !
         205             :
         206        1536 : subroutine vd_register()
         207             :
         208             :   !------------------------------------------------ !
         209             :   ! Register physics buffer fields and constituents !
         210             :   !------------------------------------------------ !
         211             :
         212        1536 :   use physics_buffer,      only : pbuf_add_field, dtype_r8, dtype_i4
         213             :   use trb_mtn_stress_cam,  only : trb_mtn_stress_register
         214             :   use beljaars_drag_cam,   only : beljaars_drag_register
         215             :   use eddy_diff_cam,       only : eddy_diff_register
         216             :
         217             :   ! Add fields to physics buffer
         218             :
         219             :   ! kvt is used by gw_drag.  only needs physpkg scope.
         220        1536 :   call pbuf_add_field('kvt', 'physpkg', dtype_r8, (/pcols,pverp/), kvt_idx)
         221             :
         222             :
         223        1536 :   if (eddy_scheme /= 'CLUBB_SGS') then
         224           0 :      call pbuf_add_field('kvh',      'global', dtype_r8, (/pcols, pverp/), kvh_idx)
         225             :   end if
         226             :
         227        1536 :   call pbuf_add_field('kvm',      'global', dtype_r8, (/pcols, pverp/), kvm_idx )
         228        1536 :   call pbuf_add_field('pblh',     'global', dtype_r8, (/pcols/),        pblh_idx)
         229        1536 :   call pbuf_add_field('tke',      'global', dtype_r8, (/pcols, pverp/), tke_idx)
         230             :
         231        1536 :   call pbuf_add_field('tauresx',  'global', dtype_r8, (/pcols/),        tauresx_idx)
         232        1536 :   call pbuf_add_field('tauresy',  'global', dtype_r8, (/pcols/),        tauresy_idx)
         233             :
         234        1536 :   call pbuf_add_field('tpert', 'global', dtype_r8, (/pcols/),                       tpert_idx)
         235        1536 :   call pbuf_add_field('qpert', 'global', dtype_r8, (/pcols,pcnst/),                 qpert_idx)
         236             :
         237        1536 :   if (trim(shallow_scheme) == 'UNICON') then
         238           0 :      call pbuf_add_field('qtl_flx',  'global', dtype_r8, (/pcols, pverp/), qtl_flx_idx)
         239           0 :      call pbuf_add_field('qti_flx',  'global', dtype_r8, (/pcols, pverp/), qti_flx_idx)
         240             :   end if
         241             :
         242             :   ! diag_TKE fields
         243        1536 :   if (eddy_scheme == 'diag_TKE' .or. eddy_scheme == 'SPCAM_m2005') then
         244           0 :      call eddy_diff_register()
         245             :   end if
         246             :
         247             :   ! TMS fields
         248        1536 :   call trb_mtn_stress_register()
         249             :
         250             :   ! Beljaars fields
         251        1536 :   call beljaars_drag_register()
         252             :
         253        1536 : end subroutine vd_register
         254             :
         255             : ! =============================================================================== !
         256             : !                                                                                 !
         257             : ! =============================================================================== !
         258             :
         259        1536 : subroutine vertical_diffusion_init(pbuf2d)
         260             :
         261             :   !------------------------------------------------------------------!
         262             :   ! Initialization of time independent fields for vertical diffusion !
         263             :   ! Calls initialization routines for subsidiary modules             !
         264             :   !----------------------------------------------------------------- !
         265             :
         266        1536 :   use cam_history,       only : addfld, add_default, horiz_only
         267             :   use cam_history,       only : register_vector_field
         268             :   use eddy_diff_cam,     only : eddy_diff_init
         269             :   use hb_diff,           only : init_hb_diff
         270             :   use molec_diff,        only : init_molec_diff
         271             :   use diffusion_solver,  only : init_vdiff, new_fieldlist_vdiff, vdiff_select
         272             :   use constituents,      only : cnst_get_ind, cnst_get_type_byind, cnst_name, cnst_get_molec_byind
         273             :   use spmd_utils,        only : masterproc
         274             :   use ref_pres,          only : press_lim_idx, pref_mid
         275             :   use physics_buffer,    only : pbuf_set_field, pbuf_get_index, physics_buffer_desc
         276             :   use rad_constituents,  only : rad_cnst_get_info, rad_cnst_get_mode_num_idx, &
         277             :        rad_cnst_get_mam_mmr_idx
         278             :   use trb_mtn_stress_cam,only : trb_mtn_stress_init
         279             :   use beljaars_drag_cam, only : beljaars_drag_init
         280             :   use upper_bc,          only : ubc_init
         281             :   use phys_control,      only : waccmx_is, fv_am_correction
         282             :   use ref_pres,          only : ptop_ref
         283             :
         284             :   type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         285             :   character(128) :: errstring   ! Error status for init_vdiff
         286             :   integer        :: ntop_eddy   ! Top    interface level to which eddy vertical diffusion is applied ( = 1 )
         287             :   integer        :: nbot_eddy   ! Bottom interface level to which eddy vertical diffusion is applied ( = pver )
         288             :   integer        :: k           ! Vertical loop index
         289             :
         290             :   real(r8), parameter :: ntop_eddy_pres = 1.e-7_r8 ! Pressure below which eddy diffusion is not done in WACCM-X. (Pa)
         291             :
         292             :   integer :: im, l, m, nmodes, nspec, ierr
         293             :
         294             :   logical :: history_amwg                 ! output the variables used by the AMWG diag package
         295             :   logical :: history_eddy                 ! output the eddy variables
         296             :   logical :: history_budget               ! Output tendencies and state variables for CAM4 T, qv, ql, qi
         297             :   integer :: history_budget_histfile_num  ! output history file number for budget fields
         298             :   logical :: history_waccm                ! output variables of interest for WACCM runs
         299             :
         300             :   !
         301             :   ! add sponge layer vertical diffusion
         302             :   !
         303        1536 :   if (ptop_ref>1e-1_r8.and.ptop_ref<100.0_r8) then
         304             :      !
         305             :      ! CAM7 FMT (but not CAM6 top (~225 Pa) or CAM7 low top or lower)
         306             :      !
         307        1536 :      allocate(kvm_sponge(4), stat=ierr)
         308        1536 :      if( ierr /= 0 ) then
         309           0 :         write(iulog,*) 'vertical_diffusion_init:  kvm_sponge allocation error = ',ierr
         310           0 :         call endrun('vertical_diffusion_init: failed to allocate kvm_sponge array')
         311             :      end if
         312        1536 :      kvm_sponge(1) = 2E6_r8
         313        1536 :      kvm_sponge(2) = 2E6_r8
         314        1536 :      kvm_sponge(3) = 0.5E6_r8
         315        1536 :      kvm_sponge(4) = 0.1E6_r8
         316           0 :   else if (ptop_ref>1e-4_r8) then
         317             :      !
         318             :      ! WACCM and WACCM-x
         319             :      !
         320           0 :      allocate(kvm_sponge(6), stat=ierr)
         321           0 :      if( ierr /= 0 ) then
         322           0 :         write(iulog,*) 'vertical_diffusion_init:  kvm_sponge allocation error = ',ierr
         323           0 :         call endrun('vertical_diffusion_init: failed to allocate kvm_sponge array')
         324             :      end if
         325           0 :      kvm_sponge(1) = 2E6_r8
         326           0 :      kvm_sponge(2) = 2E6_r8
         327           0 :      kvm_sponge(3) = 1.5E6_r8
         328           0 :      kvm_sponge(4) = 1.0E6_r8
         329           0 :      kvm_sponge(5) = 0.5E6_r8
         330           0 :      kvm_sponge(6) = 0.1E6_r8
         331             :   end if
         332             :
         333        1536 :   if (masterproc) then
         334           2 :      write(iulog,*)'Initializing vertical diffusion (vertical_diffusion_init)'
         335           2 :      if (allocated(kvm_sponge)) then
         336           2 :         write(iulog,*)'Artificial sponge layer vertical diffusion added:'
         337          10 :         do k=1,size(kvm_sponge(:),1)
         338           8 :            write(iulog,'(a44,i2,a17,e7.2,a8)') 'vertical diffusion coefficient at interface',k,' is increased by ', &
         339          18 :                                                 kvm_sponge(k),' m2 s-2'
         340             :         end do
         341             :      end if !allocated
         342             :   end if
         343             :
         344             :   ! Check to see if WACCM-X is on (currently we don't care whether the
         345             :   ! ionosphere is on or not, since this neutral diffusion code is the
         346             :   ! same either way).
         347        1536 :   waccmx_mode = waccmx_is('ionosphere') .or. waccmx_is('neutral')
         348             :
         349             :   ! ----------------------------------------------------------------- !
         350             :   ! Get indices of cloud liquid and ice within the constituents array !
         351             :   ! ----------------------------------------------------------------- !
         352             :
         353        1536 :   call cnst_get_ind( 'CLDLIQ', ixcldliq )
         354        1536 :   call cnst_get_ind( 'CLDICE', ixcldice )
         355             :   ! These are optional; with the CAM4 microphysics, there are no number
         356             :   ! constituents.
         357        1536 :   call cnst_get_ind( 'NUMLIQ', ixnumliq, abort=.false. )
         358        1536 :   call cnst_get_ind( 'NUMICE', ixnumice, abort=.false. )
         359             :
         360             :   ! prog_modal_aero determines whether prognostic modal aerosols are present in the run.
         361        1536 :   call phys_getopts(prog_modal_aero_out=prog_modal_aero)
         362        1536 :   if (prog_modal_aero) then
         363             :
         364             :      ! Get the constituent indices of the number and mass mixing ratios of the modal
         365             :      ! aerosols.
         366             :      !
         367             :      ! N.B. - This implementation assumes that the prognostic modal aerosols are
         368             :      !        impacting the climate calculation (i.e., can get info from list 0).
         369             :      !
         370             :
         371             :      ! First need total number of mam constituents
         372           0 :      call rad_cnst_get_info(0, nmodes=nmodes)
         373           0 :      do m = 1, nmodes
         374           0 :         call rad_cnst_get_info(0, m, nspec=nspec)
         375           0 :         pmam_ncnst = pmam_ncnst + 1 + nspec
         376             :      end do
         377             :
         378           0 :      allocate(pmam_cnst_idx(pmam_ncnst))
         379             :
         380             :      ! Get the constituent indicies
         381           0 :      im = 1
         382           0 :      do m = 1, nmodes
         383           0 :         call rad_cnst_get_mode_num_idx(m, pmam_cnst_idx(im))
         384           0 :         im = im + 1
         385           0 :         call rad_cnst_get_info(0, m, nspec=nspec)
         386           0 :         do l = 1, nspec
         387           0 :            call rad_cnst_get_mam_mmr_idx(m, l, pmam_cnst_idx(im))
         388           0 :            im = im + 1
         389             :         end do
         390             :      end do
         391             :   end if
         392             :
         393             :   ! Initialize upper boundary condition module
         394             :
         395        1536 :   call ubc_init()
         396             :
         397             :   ! ---------------------------------------------------------------------------------------- !
         398             :   ! Initialize molecular diffusivity module                                                  !
         399             :   ! Note that computing molecular diffusivities is a trivial expense, but constituent        !
         400             :   ! diffusivities depend on their molecular weights. Decomposing the diffusion matrix        !
         401             :   ! for each constituent is a needless expense unless the diffusivity is significant.        !
         402             :   ! ---------------------------------------------------------------------------------------- !
         403             :
         404             :   !----------------------------------------------------------------------------------------
         405             :   ! Initialize molecular diffusion and get top and bottom molecular diffusion limits
         406             :   !----------------------------------------------------------------------------------------
         407             :
         408        1536 :   if( do_molec_diff ) then
         409             :      call init_molec_diff( r8, pcnst, mwdry, avogad, &
         410           0 :           errstring)
         411             :
         412           0 :      call handle_errmsg(errstring, subname="init_molec_diff")
         413             :
         414           0 :      call addfld( 'TTPXMLC', horiz_only, 'A', 'K/S', 'Top interf. temp. flux: molec. viscosity' )
         415           0 :      if( masterproc ) write(iulog,fmt='(a,i3,5x,a,i3)') 'NBOT_MOLEC =', nbot_molec
         416             :   end if
         417             :
         418             :   ! ---------------------------------- !
         419             :   ! Initialize eddy diffusivity module !
         420             :   ! ---------------------------------- !
         421             :
         422             :   ! ntop_eddy must be 1 or <= nbot_molec
         423             :   ! Currently, it is always 1 except for WACCM-X.
         424        1536 :   if ( waccmx_mode ) then
         425           0 :      ntop_eddy  = press_lim_idx(ntop_eddy_pres, top=.true.)
         426             :   else
         427        1536 :      ntop_eddy = 1
         428             :   end if
         429        1536 :   nbot_eddy  = pver
         430             :
         431        1536 :   if (masterproc) write(iulog, fmt='(a,i3,5x,a,i3)') 'NTOP_EDDY  =', ntop_eddy, 'NBOT_EDDY  =', nbot_eddy
         432             :
         433        1536 :   call phys_getopts(do_hb_above_clubb_out=do_hb_above_clubb)
         434             :
         435           0 :   select case ( eddy_scheme )
         436             :   case ( 'diag_TKE', 'SPCAM_m2005' )
         437           0 :      if( masterproc ) write(iulog,*) &
         438           0 :           'vertical_diffusion_init: eddy_diffusivity scheme: UW Moist Turbulence Scheme by Bretherton and Park'
         439           0 :      call eddy_diff_init(pbuf2d, ntop_eddy, nbot_eddy)
         440             :   case ( 'HB', 'HBR', 'SPCAM_sam1mom')
         441           0 :      if( masterproc ) write(iulog,*) 'vertical_diffusion_init: eddy_diffusivity scheme:  Holtslag and Boville'
         442             :      call init_hb_diff(gravit, cpair, ntop_eddy, nbot_eddy, pref_mid, &
         443           0 :           karman, eddy_scheme)
         444           0 :      call addfld('HB_ri',      (/ 'lev' /),  'A', 'no',  'Richardson Number (HB Scheme), I' )
         445             :   case ( 'CLUBB_SGS' )
         446        1536 :      do_pbl_diags = .true.
         447        1536 :      call init_hb_diff(gravit, cpair, ntop_eddy, nbot_eddy, pref_mid, karman, eddy_scheme)
         448             :      !
         449             :      ! run HB scheme where CLUBB is not active when running cam7 or cam6 physics
         450             :      ! else init_hb_diff is called just for diagnostic purposes
         451             :      !
         452        3072 :      if (do_hb_above_clubb) then
         453        1536 :        if( masterproc ) then
         454           2 :          write(iulog,*) 'vertical_diffusion_init: '
         455           2 :          write(iulog,*) 'eddy_diffusivity scheme where CLUBB is not active:  Holtslag and Boville'
         456             :        end if
         457        3072 :        call addfld('HB_ri',      (/ 'lev' /),  'A', 'no',  'Richardson Number (HB Scheme), I' )
         458             :      end if
         459             :   end select
         460             :
         461             :   ! ------------------------------------------- !
         462             :   ! Initialize turbulent mountain stress module !
         463             :   ! ------------------------------------------- !
         464             :
         465        1536 :   call trb_mtn_stress_init()
         466             :
         467             :   ! ----------------------------------- !
         468             :   ! Initialize Beljaars SGO drag module !
         469             :   ! ----------------------------------- !
         470             :
         471        1536 :   call beljaars_drag_init()
         472             :
         473             :   ! ---------------------------------- !
         474             :   ! Initialize diffusion solver module !
         475             :   ! ---------------------------------- !
         476             :
         477        1536 :   call init_vdiff(r8, iulog, rair, cpair, gravit, do_iss, fv_am_correction, errstring)
         478        1536 :   call handle_errmsg(errstring, subname="init_vdiff")
         479             :
         480             :   ! Use fieldlist_wet to select the fields which will be diffused using moist mixing ratios ( all by default )
         481             :   ! Use fieldlist_dry to select the fields which will be diffused using dry   mixing ratios.
         482             :
         483        1536 :   fieldlist_wet = new_fieldlist_vdiff( pcnst)
         484        1536 :   fieldlist_dry = new_fieldlist_vdiff( pcnst)
         485        1536 :   fieldlist_molec = new_fieldlist_vdiff( pcnst)
         486             :
         487        1536 :   if( vdiff_select( fieldlist_wet, 'u' ) .ne. '' ) call endrun( vdiff_select( fieldlist_wet, 'u' ) )
         488        1536 :   if( vdiff_select( fieldlist_wet, 'v' ) .ne. '' ) call endrun( vdiff_select( fieldlist_wet, 'v' ) )
         489        1536 :   if( vdiff_select( fieldlist_wet, 's' ) .ne. '' ) call endrun( vdiff_select( fieldlist_wet, 's' ) )
         490             :
         491       18432 :   constit_loop: do k = 1, pcnst
         492             :
         493       16896 :      if (prog_modal_aero) then
         494             :         ! Do not diffuse droplet number - treated in dropmixnuc
         495           0 :         if (k == ixnumliq) cycle constit_loop
         496             :         ! Don't diffuse modal aerosol - treated in dropmixnuc
         497           0 :         do m = 1, pmam_ncnst
         498           0 :            if (k == pmam_cnst_idx(m)) cycle constit_loop
         499             :         enddo
         500             :      end if
         501             :
         502             :      ! Convert all constituents to wet before doing diffusion.
         503       16896 :      if( vdiff_select( fieldlist_wet, 'q', k ) .ne. '' ) call endrun( vdiff_select( fieldlist_wet, 'q', k ) )
         504             :
         505             :      ! ----------------------------------------------- !
         506             :      ! Select constituents for molecular diffusion     !
         507             :      ! ----------------------------------------------- !
         508       18432 :      if ( cnst_get_molec_byind(k) .eq. 'minor' ) then
         509       33792 :         if( vdiff_select(fieldlist_molec,'q',k) .ne. '' ) call endrun( vdiff_select( fieldlist_molec,'q',k ) )
         510             :      endif
         511             :
         512             :   end do constit_loop
         513             :
         514             :   ! ------------------------ !
         515             :   ! Diagnostic output fields !
         516             :   ! ------------------------ !
         517             :
         518       18432 :   do k = 1, pcnst
         519       16896 :      vdiffnam(k) = 'VD'//cnst_name(k)
         520       16896 :      if( k == 1 ) vdiffnam(k) = 'VD01'    !**** compatibility with old code ****
         521       35328 :      call addfld( vdiffnam(k), (/ 'lev' /), 'A', 'kg/kg/s', 'Vertical diffusion of '//cnst_name(k) )
         522             :   end do
         523             :
         524        1536 :   if (.not. do_pbl_diags) then
         525           0 :      call addfld( 'PBLH'        , horiz_only    , 'A', 'm'      , 'PBL height'                                         )
         526           0 :      call addfld( 'QT'          , (/ 'lev' /)   , 'A', 'kg/kg'  , 'Total water mixing ratio'                           )
         527           0 :      call addfld( 'SL'          , (/ 'lev' /)   , 'A', 'J/kg'   , 'Liquid water static energy'                         )
         528           0 :      call addfld( 'SLV'         , (/ 'lev' /)   , 'A', 'J/kg'   , 'Liq wat virtual static energy'                      )
         529           0 :      call addfld( 'SLFLX'       , (/ 'ilev' /)  , 'A', 'W/m2'   , 'Liquid static energy flux'                          )
         530           0 :      call addfld( 'QTFLX'       , (/ 'ilev' /)  , 'A', 'W/m2'   , 'Total water flux'                                   )
         531           0 :      call addfld( 'TKE'         , (/ 'ilev' /)  , 'A', 'm2/s2'  , 'Turbulent Kinetic Energy'                           )
         532           0 :      call addfld( 'TPERT'       , horiz_only    , 'A', 'K'      , 'Perturbation temperature (eddies in PBL)'           )
         533           0 :      call addfld( 'QPERT'       , horiz_only    , 'A', 'kg/kg'  , 'Perturbation specific humidity (eddies in PBL)'     )
         534             :
         535           0 :      call addfld( 'UFLX'        , (/ 'ilev' /)  , 'A', 'W/m2'   , 'Zonal momentum flux'                                )
         536           0 :      call addfld( 'VFLX'        , (/ 'ilev' /)  , 'A', 'W/m2'   , 'Meridional momentm flux'                            )
         537           0 :      call register_vector_field('UFLX', 'VFLX')
         538             :   end if
         539             :
         540        1536 :   call addfld( 'USTAR'       , horiz_only    , 'A', 'm/s'    , 'Surface friction velocity'                          )
         541        3072 :   call addfld( 'KVH'         , (/ 'ilev' /)  , 'A', 'm2/s'   , 'Vertical diffusion diffusivities (heat/moisture)'   )
         542        3072 :   call addfld( 'KVM'         , (/ 'ilev' /)  , 'A', 'm2/s'   , 'Vertical diffusion diffusivities (momentum)'        )
         543        3072 :   call addfld( 'KVT'         , (/ 'ilev' /)  , 'A', 'm2/s'   , 'Vertical diffusion kinematic molecular conductivity')
         544        3072 :   call addfld( 'CGS'         , (/ 'ilev' /)  , 'A', 's/m2'   , 'Counter-gradient coeff on surface kinematic fluxes' )
         545        3072 :   call addfld( 'DTVKE'       , (/ 'lev' /)   , 'A', 'K/s'    , 'dT/dt vertical diffusion KE dissipation'            )
         546        3072 :   call addfld( 'DTV'         , (/ 'lev' /)   , 'A', 'K/s'    , 'T vertical diffusion'                               )
         547        3072 :   call addfld( 'DUV'         , (/ 'lev' /)   , 'A', 'm/s2'   , 'U vertical diffusion'                               )
         548        3072 :   call addfld( 'DVV'         , (/ 'lev' /)   , 'A', 'm/s2'   , 'V vertical diffusion'                               )
         549             :
         550             :   ! ---------------------------------------------------------------------------- !
         551             :   ! Below ( with '_PBL') are for detailed analysis of UW Moist Turbulence Scheme !
         552             :   ! ---------------------------------------------------------------------------- !
         553             :
         554        1536 :   if (.not. do_pbl_diags) then
         555             :
         556           0 :      call addfld( 'qt_pre_PBL',   (/ 'lev' /)   , 'A', 'kg/kg'  , 'qt_prePBL'          )
         557           0 :      call addfld( 'sl_pre_PBL',   (/ 'lev' /)   , 'A', 'J/kg'   , 'sl_prePBL'          )
         558           0 :      call addfld( 'slv_pre_PBL',  (/ 'lev' /)   , 'A', 'J/kg'   , 'slv_prePBL'         )
         559           0 :      call addfld( 'u_pre_PBL',    (/ 'lev' /)   , 'A', 'm/s'    , 'u_prePBL'           )
         560           0 :      call addfld( 'v_pre_PBL',    (/ 'lev' /)   , 'A', 'm/s'    , 'v_prePBL'           )
         561           0 :      call addfld( 'qv_pre_PBL',   (/ 'lev' /)   , 'A', 'kg/kg'  , 'qv_prePBL'          )
         562           0 :      call addfld( 'ql_pre_PBL',   (/ 'lev' /)   , 'A', 'kg/kg'  , 'ql_prePBL'          )
         563           0 :      call addfld( 'qi_pre_PBL',   (/ 'lev' /)   , 'A', 'kg/kg'  , 'qi_prePBL'          )
         564           0 :      call addfld( 't_pre_PBL',    (/ 'lev' /)   , 'A', 'K'      , 't_prePBL'           )
         565           0 :      call addfld( 'rh_pre_PBL',   (/ 'lev' /)   , 'A', '%'      , 'rh_prePBL'          )
         566             :
         567           0 :      call addfld( 'qt_aft_PBL',   (/ 'lev' /)   , 'A', 'kg/kg'  , 'qt_afterPBL'        )
         568           0 :      call addfld( 'sl_aft_PBL',   (/ 'lev' /)   , 'A', 'J/kg'   , 'sl_afterPBL'        )
         569           0 :      call addfld( 'slv_aft_PBL',  (/ 'lev' /)   , 'A', 'J/kg'   , 'slv_afterPBL'       )
         570           0 :      call addfld( 'u_aft_PBL',    (/ 'lev' /)   , 'A', 'm/s'    , 'u_afterPBL'         )
         571           0 :      call addfld( 'v_aft_PBL',    (/ 'lev' /)   , 'A', 'm/s'    , 'v_afterPBL'         )
         572           0 :      call addfld( 'qv_aft_PBL',   (/ 'lev' /)   , 'A', 'kg/kg'  , 'qv_afterPBL'        )
         573           0 :      call addfld( 'ql_aft_PBL',   (/ 'lev' /)   , 'A', 'kg/kg'  , 'ql_afterPBL'        )
         574           0 :      call addfld( 'qi_aft_PBL',   (/ 'lev' /)   , 'A', 'kg/kg'  , 'qi_afterPBL'        )
         575           0 :      call addfld( 't_aft_PBL',    (/ 'lev' /)   , 'A', 'K'      , 't_afterPBL'         )
         576           0 :      call addfld( 'rh_aft_PBL',   (/ 'lev' /)   , 'A', '%'      , 'rh_afterPBL'        )
         577             :
         578           0 :      call addfld( 'slflx_PBL',    (/ 'ilev' /)  , 'A', 'J/m2/s' , 'sl flux by PBL'     )
         579           0 :      call addfld( 'qtflx_PBL',    (/ 'ilev' /)  , 'A', 'kg/m2/s', 'qt flux by PBL'     )
         580           0 :      call addfld( 'uflx_PBL',     (/ 'ilev' /)  , 'A', 'kg/m/s2', 'u flux by PBL'      )
         581           0 :      call addfld( 'vflx_PBL',     (/ 'ilev' /)  , 'A', 'kg/m/s2', 'v flux by PBL'      )
         582             :
         583           0 :      call addfld( 'slflx_cg_PBL', (/ 'ilev' /)  , 'A', 'J/m2/s' , 'sl_cg flux by PBL'  )
         584           0 :      call addfld( 'qtflx_cg_PBL', (/ 'ilev' /)  , 'A', 'kg/m2/s', 'qt_cg flux by PBL'  )
         585           0 :      call addfld( 'uflx_cg_PBL',  (/ 'ilev' /)  , 'A', 'kg/m/s2', 'u_cg flux by PBL'   )
         586           0 :      call addfld( 'vflx_cg_PBL',  (/ 'ilev' /)  , 'A', 'kg/m/s2', 'v_cg flux by PBL'   )
         587             :
         588           0 :      call addfld( 'qtten_PBL',    (/ 'lev' /)   , 'A', 'kg/kg/s', 'qt tendency by PBL' )
         589           0 :      call addfld( 'slten_PBL',    (/ 'lev' /)   , 'A', 'J/kg/s' , 'sl tendency by PBL' )
         590           0 :      call addfld( 'uten_PBL',     (/ 'lev' /)   , 'A', 'm/s2'   , 'u tendency by PBL'  )
         591           0 :      call addfld( 'vten_PBL',     (/ 'lev' /)   , 'A', 'm/s2'   , 'v tendency by PBL'  )
         592           0 :      call addfld( 'qvten_PBL',    (/ 'lev' /)   , 'A', 'kg/kg/s', 'qv tendency by PBL' )
         593           0 :      call addfld( 'qlten_PBL',    (/ 'lev' /)   , 'A', 'kg/kg/s', 'ql tendency by PBL' )
         594           0 :      call addfld( 'qiten_PBL',    (/ 'lev' /)   , 'A', 'kg/kg/s', 'qi tendency by PBL' )
         595           0 :      call addfld( 'tten_PBL',     (/ 'lev' /)   , 'A', 'K/s'    , 'T tendency by PBL'  )
         596           0 :      call addfld( 'rhten_PBL',    (/ 'lev' /)   , 'A', '%/s'    , 'RH tendency by PBL' )
         597             :
         598             :   end if
         599             :
         600        1536 :   call addfld ('ustar',horiz_only, 'A',     ' ',' ')
         601        1536 :   call addfld ('obklen',horiz_only, 'A',    ' ',' ')
         602             :
         603             :   ! ----------------------------
         604             :   ! determine default variables
         605             :   ! ----------------------------
         606             :
         607             :   call phys_getopts( history_amwg_out = history_amwg, &
         608             :        history_eddy_out = history_eddy, &
         609             :        history_budget_out = history_budget, &
         610             :        history_budget_histfile_num_out = history_budget_histfile_num, &
         611        1536 :        history_waccm_out = history_waccm)
         612             :
         613        1536 :   if (history_amwg) then
         614        1536 :      call add_default(  vdiffnam(1), 1, ' ' )
         615        1536 :      call add_default( 'DTV'       , 1, ' ' )
         616        1536 :      if (.not. do_pbl_diags) then
         617           0 :         call add_default( 'PBLH'      , 1, ' ' )
         618             :      end if
         619             :   endif
         620             :
         621        1536 :   if (history_eddy) then
         622           0 :      call add_default( 'UFLX    ', 1, ' ' )
         623           0 :      call add_default( 'VFLX    ', 1, ' ' )
         624             :   endif
         625             :
         626        1536 :   if( history_budget ) then
         627           0 :      call add_default( vdiffnam(ixcldliq), history_budget_histfile_num, ' ' )
         628           0 :      call add_default( vdiffnam(ixcldice), history_budget_histfile_num, ' ' )
         629           0 :      if( history_budget_histfile_num > 1 ) then
         630           0 :         call add_default(  vdiffnam(1), history_budget_histfile_num, ' ' )
         631           0 :         call add_default( 'DTV'       , history_budget_histfile_num, ' ' )
         632             :      end if
         633             :   end if
         634             :
         635        1536 :   if ( history_waccm ) then
         636           0 :      if (do_molec_diff) then
         637           0 :         call add_default ( 'TTPXMLC', 1, ' ' )
         638             :      end if
         639           0 :      call add_default( 'DUV'     , 1, ' ' )
         640           0 :      call add_default( 'DVV'     , 1, ' ' )
         641             :   end if
         642             :   ! ----------------------------
         643             :
         644             :
         645        1536 :   ksrftms_idx = pbuf_get_index('ksrftms')
         646        1536 :   tautmsx_idx = pbuf_get_index('tautmsx')
         647        1536 :   tautmsy_idx = pbuf_get_index('tautmsy')
         648             :
         649        1536 :   dragblj_idx = pbuf_get_index('dragblj')
         650        1536 :   taubljx_idx = pbuf_get_index('taubljx')
         651        1536 :   taubljy_idx = pbuf_get_index('taubljy')
         652             :
         653        1536 :   if (eddy_scheme == 'CLUBB_SGS') then
         654        1536 :      kvh_idx = pbuf_get_index('kvh')
         655             :   end if
         656             :
         657        1536 :   if (do_hb_above_clubb) then
         658             :      ! pbuf field denoting top of clubb
         659        1536 :      clubbtop_idx = pbuf_get_index('clubbtop')
         660             :   end if
         661             :
         662             :   ! Initialization of some pbuf fields
         663        1536 :   if (is_first_step()) then
         664             :      ! Initialization of pbuf fields tke, kvh, kvm are done in phys_inidat
         665         768 :      call pbuf_set_field(pbuf2d, tauresx_idx,  0.0_r8)
         666         768 :      call pbuf_set_field(pbuf2d, tauresy_idx,  0.0_r8)
         667         768 :      if (trim(shallow_scheme) == 'UNICON') then
         668           0 :         call pbuf_set_field(pbuf2d, qtl_flx_idx,  0.0_r8)
         669           0 :         call pbuf_set_field(pbuf2d, qti_flx_idx,  0.0_r8)
         670             :      end if
         671             :   end if
         672        1536 : end subroutine vertical_diffusion_init
         673             :
         674             : ! =============================================================================== !
         675             : !                                                                                 !
         676             : ! =============================================================================== !
         677             :
         678      741888 : subroutine vertical_diffusion_ts_init( pbuf2d, state )
         679             :
         680             :   !-------------------------------------------------------------- !
         681             :   ! Timestep dependent setting,                                   !
         682             :   ! At present only invokes upper bc code                         !
         683             :   !-------------------------------------------------------------- !
         684        1536 :   use upper_bc,       only : ubc_timestep_init
         685             :   use physics_types , only : physics_state
         686             :   use ppgrid        , only : begchunk, endchunk
         687             :
         688             :   use physics_buffer, only : physics_buffer_desc
         689             :
         690             :   type(physics_state), intent(in) :: state(begchunk:endchunk)
         691             :   type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         692             :
         693      370944 :   call ubc_timestep_init( pbuf2d, state)
         694             :
         695      370944 : end subroutine vertical_diffusion_ts_init
         696             :
         697             : ! =============================================================================== !
         698             : !                                                                                 !
         699             : ! =============================================================================== !
         700             :
         701    17870112 : subroutine vertical_diffusion_tend( &
         702             :      ztodt    , state    , cam_in,          &
         703             :      ustar    , obklen   , ptend    , &
         704             :      cldn     , pbuf)
         705             :   !---------------------------------------------------- !
         706             :   ! This is an interface routine for vertical diffusion !
         707             :   !---------------------------------------------------- !
         708      370944 :   use physics_buffer,     only : physics_buffer_desc, pbuf_get_field, pbuf_set_field
         709             :   use physics_types,      only : physics_state, physics_ptend, physics_ptend_init
         710             :   use physics_types,      only : set_dry_to_wet, set_wet_to_dry
         711             :
         712             :   use camsrfexch,         only : cam_in_t
         713             :   use cam_history,        only : outfld
         714             :
         715             :   use trb_mtn_stress_cam, only : trb_mtn_stress_tend
         716             :   use beljaars_drag_cam,  only : beljaars_drag_tend
         717             :   use eddy_diff_cam,      only : eddy_diff_tend
         718             :   use hb_diff,            only : compute_hb_diff, compute_hb_free_atm_diff
         719             :   use wv_saturation,      only : qsat
         720             :   use molec_diff,         only : compute_molec_diff, vd_lu_qdecomp
         721             :   use constituents,       only : qmincg, qmin, cnst_type
         722             :   use diffusion_solver,   only : compute_vdiff, any, operator(.not.)
         723             :   use air_composition,    only : cpairv, rairv !Needed for calculation of upward H flux
         724             :   use time_manager,       only : get_nstep
         725             :   use constituents,       only : cnst_get_type_byind, cnst_name, &
         726             :        cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx
         727             :   use physconst,          only : pi
         728             :   use pbl_utils,          only : virtem, calc_obklen, calc_ustar
         729             :   use upper_bc,           only : ubc_get_vals, ubc_fixed_temp
         730             :   use upper_bc,           only : ubc_get_flxs
         731             :   use coords_1d,          only : Coords1D
         732             :   use phys_control,       only : cam_physpkg_is
         733             :   use ref_pres,           only : ptop_ref
         734             :
         735             :   ! --------------- !
         736             :   ! Input Arguments !
         737             :   ! --------------- !
         738             :
         739             :   type(physics_state), intent(inout) :: state                     ! Physics state variables
         740             :   type(cam_in_t),      intent(in)    :: cam_in                    ! Surface inputs
         741             :
         742             :   real(r8),            intent(in)    :: ztodt                     ! 2 delta-t [ s ]
         743             :   real(r8),            intent(in)    :: cldn(pcols,pver)          ! New stratus fraction [ fraction ]
         744             :
         745             :   ! ---------------------- !
         746             :   ! Input-Output Arguments !
         747             :   ! ---------------------- !
         748             :
         749             :   type(physics_ptend), intent(out) :: ptend                       ! Individual parameterization tendencies
         750             :   type(physics_buffer_desc), pointer :: pbuf(:)
         751             :
         752             :   ! ---------------- !
         753             :   ! Output Arguments !
         754             :   ! ---------------- !
         755             :
         756             :   real(r8),            intent(out)   :: ustar(pcols)              ! Surface friction velocity [ m/s ]
         757             :   real(r8),            intent(out)   :: obklen(pcols)             ! Obukhov length [ m ]
         758             :
         759             :   ! --------------- !
         760             :   ! Local Variables !
         761             :   ! --------------- !
         762             :
         763             :   character(128) :: errstring                                     ! Error status for compute_vdiff
         764             :
         765             :   integer  :: lchnk                                               ! Chunk identifier
         766             :   integer  :: ncol                                                ! Number of atmospheric columns
         767             :   integer  :: i, k, l, m                                          ! column, level, constituent indices
         768             :
         769             :   real(r8) :: dtk(pcols,pver)                                     ! T tendency from KE dissipation
         770     1489176 :   real(r8), pointer   :: tke(:,:)                                 ! Turbulent kinetic energy [ m2/s2 ]
         771             :
         772     1489176 :   real(r8), pointer   :: qtl_flx(:,:)                             ! overbar(w'qtl') where qtl = qv + ql
         773     1489176 :   real(r8), pointer   :: qti_flx(:,:)                             ! overbar(w'qti') where qti = qv + qi
         774             :
         775             :   real(r8) :: cgs(pcols,pverp)                                    ! Counter-gradient star  [ cg/flux ]
         776             :   real(r8) :: cgh(pcols,pverp)                                    ! Counter-gradient term for heat
         777             :   real(r8) :: rztodt                                              ! 1./ztodt [ 1/s ]
         778     1489176 :   real(r8), pointer :: ksrftms(:)                                 ! Turbulent mountain stress surface drag coefficient [ kg/s/m2 ]
         779     1489176 :   real(r8), pointer :: tautmsx(:)                                 ! U component of turbulent mountain stress [ N/m2 ]
         780     1489176 :   real(r8), pointer :: tautmsy(:)                                 ! V component of turbulent mountain stress [ N/m2 ]
         781             :   real(r8) :: tautotx(pcols)                                      ! U component of total surface stress [ N/m2 ]
         782             :   real(r8) :: tautoty(pcols)                                      ! V component of total surface stress [ N/m2 ]
         783             :
         784     1489176 :   real(r8), pointer :: dragblj(:,:)                               ! Beljaars SGO form drag profile [ 1/s ]
         785     1489176 :   real(r8), pointer :: taubljx(:)                                 ! U component of turbulent mountain stress [ N/m2 ]
         786     1489176 :   real(r8), pointer :: taubljy(:)                                 ! V component of turbulent mountain stress [ N/m2 ]
         787             :
         788     1489176 :   real(r8), pointer :: kvh_in(:,:)                                ! kvh from previous timestep [ m2/s ]
         789     1489176 :   real(r8), pointer :: kvm_in(:,:)                                ! kvm from previous timestep [ m2/s ]
         790     1489176 :   real(r8), pointer :: kvt(:,:)                                   ! Molecular kinematic conductivity for temperature [  ]
         791             :   real(r8) :: kvq(pcols,pverp)                                    ! Eddy diffusivity for constituents [ m2/s ]
         792             :   real(r8) :: kvh(pcols,pverp)                                    ! Eddy diffusivity for heat [ m2/s ]
         793             :   real(r8) :: kvm(pcols,pverp)                                    ! Eddy diffusivity for momentum [ m2/s ]
         794             :   real(r8) :: kvm_temp(pcols,pverp)                               ! Dummy eddy diffusivity for momentum (unused) [ m2/s ]
         795             :   real(r8) :: dtk_temp(pcols,pverp)                               ! Unused output from second compute_vdiff call
         796             :   real(r8) :: tautmsx_temp(pcols)                                 ! Unused output from second compute_vdiff call
         797             :   real(r8) :: tautmsy_temp(pcols)                                 ! Unused output from second compute_vdiff call
         798             :   real(r8) :: topflx_temp(pcols)                                  ! Unused output from second compute_vdiff call
         799             :   real(r8) :: sprod(pcols,pverp)                                  ! Shear production of tke [ m2/s3 ]
         800             :   real(r8) :: sfi(pcols,pverp)                                    ! Saturation fraction at interfaces [ fraction ]
         801             :   real(r8) :: sl(pcols,pver)
         802             :   real(r8) :: qt(pcols,pver)
         803             :   real(r8) :: slv(pcols,pver)
         804             :   real(r8) :: sl_prePBL(pcols,pver)
         805             :   real(r8) :: qt_prePBL(pcols,pver)
         806             :   real(r8) :: slv_prePBL(pcols,pver)
         807             :   real(r8) :: slten(pcols,pver)
         808             :   real(r8) :: qtten(pcols,pver)
         809             :   real(r8) :: slflx(pcols,pverp)
         810             :   real(r8) :: qtflx(pcols,pverp)
         811             :   real(r8) :: uflx(pcols,pverp)
         812             :   real(r8) :: vflx(pcols,pverp)
         813             :   real(r8) :: slflx_cg(pcols,pverp)
         814             :   real(r8) :: qtflx_cg(pcols,pverp)
         815             :   real(r8) :: uflx_cg(pcols,pverp)
         816             :   real(r8) :: vflx_cg(pcols,pverp)
         817             :   real(r8) :: th(pcols,pver)                                      ! Potential temperature
         818             :   real(r8) :: topflx(pcols)                                       ! Molecular heat flux at top interface
         819             :   real(r8) :: rhoair
         820             :
         821             :   real(r8) :: ri(pcols,pver)                                      ! richardson number (HB output)
         822             :
         823             :   ! for obklen calculation outside HB
         824             :   real(r8) :: thvs(pcols)                                         ! Virtual potential temperature at surface
         825             :   real(r8) :: rrho(pcols)                                         ! Reciprocal of density at surface
         826             :   real(r8) :: khfs(pcols)                                         ! sfc kinematic heat flux [mK/s]
         827             :   real(r8) :: kqfs(pcols)                                         ! sfc kinematic water vapor flux [m/s]
         828             :   real(r8) :: kbfs(pcols)                                         ! sfc kinematic buoyancy flux [m^2/s^3]
         829             :
         830             :   real(r8) :: ftem(pcols,pver)                                    ! Saturation vapor pressure before PBL
         831             :   real(r8) :: ftem_prePBL(pcols,pver)                             ! Saturation vapor pressure before PBL
         832             :   real(r8) :: ftem_aftPBL(pcols,pver)                             ! Saturation vapor pressure after PBL
         833             :   real(r8) :: tem2(pcols,pver)                                    ! Saturation specific humidity and RH
         834             :   real(r8) :: t_aftPBL(pcols,pver)                                ! Temperature after PBL diffusion
         835             :   real(r8) :: tten(pcols,pver)                                    ! Temperature tendency by PBL diffusion
         836             :   real(r8) :: rhten(pcols,pver)                                   ! RH tendency by PBL diffusion
         837             :   real(r8) :: qv_aft_PBL(pcols,pver)                              ! qv after PBL diffusion
         838             :   real(r8) :: ql_aft_PBL(pcols,pver)                              ! ql after PBL diffusion
         839             :   real(r8) :: qi_aft_PBL(pcols,pver)                              ! qi after PBL diffusion
         840             :   real(r8) :: s_aft_PBL(pcols,pver)                               ! s after PBL diffusion
         841             :   real(r8) :: u_aft_PBL(pcols,pver)                               ! u after PBL diffusion
         842             :   real(r8) :: v_aft_PBL(pcols,pver)                               ! v after PBL diffusion
         843             :   real(r8) :: qv_pro(pcols,pver)
         844             :   real(r8) :: ql_pro(pcols,pver)
         845             :   real(r8) :: qi_pro(pcols,pver)
         846             :   real(r8) :: s_pro(pcols,pver)
         847             :   real(r8) :: t_pro(pcols,pver)
         848     1489176 :   real(r8), pointer :: tauresx(:)                                      ! Residual stress to be added in vdiff to correct
         849     1489176 :   real(r8), pointer :: tauresy(:)                                      ! for turb stress mismatch between sfc and atm accumulated.
         850             :
         851             :   ! Interpolated interface values.
         852             :   real(r8) :: tint(pcols,pver+1)      ! Temperature [ K ]
         853             :   real(r8) :: rairi(pcols,pver+1)     ! Gas constant [ J/K/kg ]
         854             :   real(r8) :: rhoi(pcols,pver+1)      ! Density of air [ kg/m^3 ]
         855             :   real(r8) :: rhoi_dry(pcols,pver+1)  ! Density of air based on dry air pressure [ kg/m^3 ]
         856             :
         857             :   ! Upper boundary conditions
         858             :   real(r8) :: ubc_t(pcols)            ! Temperature [ K ]
         859             :   real(r8) :: ubc_mmr(pcols,pcnst)    ! Mixing ratios [ kg/kg ]
         860             :   real(r8) :: ubc_flux(pcols,pcnst)   ! Constituent upper boundary flux (kg/s/m^2)
         861             :
         862             :   ! Pressure coordinates used by the solver.
         863     1489176 :   type(Coords1D) :: p
         864     1489176 :   type(Coords1D) :: p_dry
         865             :
         866     1489176 :   real(r8), pointer :: tpert(:)
         867     1489176 :   real(r8), pointer :: qpert(:)
         868     1489176 :   real(r8), pointer :: pblh(:)
         869             :
         870             :   real(r8) :: tmp1(pcols)                                         ! Temporary storage
         871             :
         872             :   integer  :: nstep
         873             :   real(r8) :: sum1, sum2, sum3, pdelx
         874             :   real(r8) :: sflx
         875             :
         876             :   ! Copy state so we can pass to intent(inout) routines that return
         877             :   ! new state instead of a tendency.
         878             :   real(r8) :: s_tmp(pcols,pver)
         879             :   real(r8) :: u_tmp(pcols,pver)
         880             :   real(r8) :: v_tmp(pcols,pver)
         881             :   real(r8) :: q_tmp(pcols,pver,pcnst)
         882             :
         883             :   ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity
         884             :   real(r8) :: kq_scal(pcols,pver+1)
         885             :   ! composition dependent mw_fac on interface level
         886             :   real(r8) :: mw_fac(pcols,pver+1,pcnst)
         887             :
         888             :   ! Dry static energy top boundary condition.
         889             :   real(r8) :: dse_top(pcols)
         890             :
         891             :   ! Copies of flux arrays used to zero out any parts that are applied
         892             :   ! elsewhere (e.g. by CLUBB).
         893             :   real(r8) :: taux(pcols)
         894             :   real(r8) :: tauy(pcols)
         895             :   real(r8) :: shflux(pcols)
         896             :   real(r8) :: cflux(pcols,pcnst)
         897     1489176 :   integer,  pointer :: clubbtop(:)   ! (pcols)
         898             :
         899             :   logical  :: lq(pcnst)
         900             :
         901             :   ! ----------------------- !
         902             :   ! Main Computation Begins !
         903             :   ! ----------------------- !
         904             :
         905             :   ! Assume 'wet' mixing ratios in diffusion code.
         906     1489176 :   call set_dry_to_wet(state, convert_cnst_type='dry')
         907             :
         908     1489176 :   rztodt = 1._r8 / ztodt
         909     1489176 :   lchnk  = state%lchnk
         910     1489176 :   ncol   = state%ncol
         911             :
         912     1489176 :   call pbuf_get_field(pbuf, tauresx_idx,  tauresx)
         913     1489176 :   call pbuf_get_field(pbuf, tauresy_idx,  tauresy)
         914     1489176 :   call pbuf_get_field(pbuf, tpert_idx,    tpert)
         915     1489176 :   call pbuf_get_field(pbuf, qpert_idx,    qpert)
         916     1489176 :   call pbuf_get_field(pbuf, pblh_idx,     pblh)
         917             :
         918             :   ! Interpolate temperature to interfaces.
         919   138493368 :   do k = 2, pver
         920  2289140568 :      do i = 1, ncol
         921  2287651392 :         tint(i,k)  = 0.5_r8 * ( state%t(i,k) + state%t(i,k-1) )
         922             :      end do
         923             :   end do
         924    24865776 :   tint(:ncol,pver+1) = state%t(:ncol,pver)
         925             :
         926             :   ! Get upper boundary values
         927     1489176 :   call ubc_get_vals( state%lchnk, ncol, state%pint, state%zi, ubc_t, ubc_mmr )
         928             :
         929     1489176 :   if (waccmx_mode) then
         930           0 :      call ubc_get_flxs( state%lchnk, ncol, state%pint, state%zi, state%t, state%q, state%phis, ubc_flux )
         931             :      ! For WACCM-X, set ubc temperature to extrapolate from next two lower interface level temperatures
         932           0 :      tint(:ncol,1) = 1.5_r8*tint(:ncol,2)-.5_r8*tint(:ncol,3)
         933     1489176 :   else if(ubc_fixed_temp) then
         934           0 :      tint(:ncol,1) = ubc_t(:ncol)
         935             :   else
         936    24865776 :      tint(:ncol,1) = state%t(:ncol,1)
         937             :   end if
         938             :
         939             :   ! Set up pressure coordinates for solver calls.
         940     1489176 :   p = Coords1D(state%pint(:ncol,:))
         941     1489176 :   p_dry = Coords1D(state%pintdry(:ncol,:))
         942             :
         943             :   !------------------------------------------------------------------------
         944             :   !  Check to see if constituent dependent gas constant needed (WACCM-X)
         945             :   !------------------------------------------------------------------------
         946     1489176 :   if (waccmx_mode) then
         947           0 :      rairi(:ncol,1) = rairv(:ncol,1,lchnk)
         948           0 :      do k = 2, pver
         949           0 :         do i = 1, ncol
         950           0 :            rairi(i,k) = 0.5_r8 * (rairv(i,k,lchnk)+rairv(i,k-1,lchnk))
         951             :         end do
         952             :      end do
         953           0 :      rairi(:ncol,pver+1) = rairv(:ncol,pver,lchnk)
         954             :   else
         955  2338872120 :      rairi(:ncol,:pver+1) = rair
         956             :   endif
         957             :
         958             :   ! Compute rho at interfaces.
         959   141471720 :   do k = 1, pver+1
         960  2338872120 :      do i = 1, ncol
         961  2337382944 :         rhoi(i,k)  = p%ifc(i,k) / (rairi(i,k)*tint(i,k))
         962             :      end do
         963             :   end do
         964             :
         965             :   ! Compute rho_dry at interfaces.
         966   141471720 :   do k = 1, pver+1
         967  2338872120 :      do i = 1, ncol
         968  2337382944 :         rhoi_dry(i,k)  = p_dry%ifc(i,k) / (rairi(i,k)*tint(i,k))
         969             :      end do
         970             :   end do
         971             :
         972             :   ! ---------------------------------------- !
         973             :   ! Computation of turbulent mountain stress !
         974             :   ! ---------------------------------------- !
         975             :
         976             :   ! Consistent with the computation of 'normal' drag coefficient, we are using
         977             :   ! the raw input (u,v) to compute 'ksrftms', not the provisionally-marched 'u,v'
         978             :   ! within the iteration loop of the PBL scheme.
         979             :
         980     1489176 :   call trb_mtn_stress_tend(state, pbuf, cam_in)
         981             :
         982     1489176 :   call pbuf_get_field(pbuf, ksrftms_idx, ksrftms)
         983     1489176 :   call pbuf_get_field(pbuf, tautmsx_idx, tautmsx)
         984     1489176 :   call pbuf_get_field(pbuf, tautmsy_idx, tautmsy)
         985             :
         986    24865776 :   tautotx(:ncol) = cam_in%wsx(:ncol) + tautmsx(:ncol)
         987    24865776 :   tautoty(:ncol) = cam_in%wsy(:ncol) + tautmsy(:ncol)
         988             :
         989             :   ! ------------------------------------- !
         990             :   ! Computation of Beljaars SGO form drag !
         991             :   ! ------------------------------------- !
         992             :
         993     1489176 :   call beljaars_drag_tend(state, pbuf, cam_in)
         994             :
         995     1489176 :   call pbuf_get_field(pbuf, dragblj_idx, dragblj)
         996     1489176 :   call pbuf_get_field(pbuf, taubljx_idx, taubljx)
         997     1489176 :   call pbuf_get_field(pbuf, taubljy_idx, taubljy)
         998             :
         999             :   ! Add Beljaars integrated drag
        1000             :
        1001    24865776 :   tautotx(:ncol) = tautotx(:ncol) + taubljx(:ncol)
        1002    24865776 :   tautoty(:ncol) = tautoty(:ncol) + taubljy(:ncol)
        1003             :
        1004             :   !----------------------------------------------------------------------- !
        1005             :   !   Computation of eddy diffusivities - Select appropriate PBL scheme    !
        1006             :   !----------------------------------------------------------------------- !
        1007     1489176 :   call pbuf_get_field(pbuf, kvm_idx,  kvm_in)
        1008     1489176 :   call pbuf_get_field(pbuf, kvh_idx,  kvh_in)
        1009     1489176 :   call pbuf_get_field(pbuf, tke_idx,  tke)
        1010             :
        1011             :   ! Get potential temperature.
        1012  2314006344 :   th(:ncol,:pver) = state%t(:ncol,:pver) * state%exner(:ncol,:pver)
        1013             :
        1014           0 :   select case (eddy_scheme)
        1015             :   case ( 'diag_TKE', 'SPCAM_m2005' )
        1016             :
        1017             :      call eddy_diff_tend(state, pbuf, cam_in, &
        1018             :           ztodt, p, tint, rhoi, cldn, wstarent, &
        1019             :           kvm_in, kvh_in, ksrftms, dragblj, tauresx, tauresy, &
        1020             :           rrho, ustar, pblh, kvm, kvh, kvq, cgh, cgs, tpert, qpert, &
        1021           0 :           tke, sprod, sfi)
        1022             :
        1023             :      ! The diag_TKE scheme does not calculate the Monin-Obukhov length, which is used in dry deposition calculations.
        1024             :      ! Use the routines from pbl_utils to accomplish this. Assumes ustar and rrho have been set.
        1025           0 :      call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol))
        1026           0 :      call calc_obklen(ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), &
        1027             :           cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol), &
        1028           0 :           khfs(:ncol),    kqfs(:ncol), kbfs(:ncol),   obklen(:ncol))
        1029             :
        1030             :
        1031             :   case ( 'HB', 'HBR', 'SPCAM_sam1mom' )
        1032             :
        1033             :      ! Modification : We may need to use 'taux' instead of 'tautotx' here, for
        1034             :      !                consistency with the previous HB scheme.
        1035             :
        1036             :
        1037             :      call compute_hb_diff(ncol                                    , &
        1038             :           th        , state%t  , state%q , state%zm , state%zi    , &
        1039             :           state%pmid, state%u  , state%v , tautotx  , tautoty     , &
        1040             :           cam_in%shf, cam_in%cflx(:,1), obklen  , ustar    , pblh , &
        1041             :           kvm       , kvh      , kvq     , cgh      , cgs         , &
        1042             :           tpert     , qpert    , cldn    , cam_in%ocnfrac  , tke  , &
        1043             :           ri        , &
        1044           0 :           eddy_scheme)
        1045             :
        1046           0 :      call outfld( 'HB_ri',          ri,         pcols,   lchnk )
        1047             :
        1048             :   case ( 'CLUBB_SGS' )
        1049             :     !
        1050             :     ! run HB scheme where CLUBB is not active when running cam7
        1051             :     !
        1052     1489176 :     if (do_hb_above_clubb) then
        1053             :       call compute_hb_free_atm_diff( ncol          , &
        1054             :            th        , state%t  , state%q , state%zm           , &
        1055             :            state%pmid, state%u  , state%v , tautotx  , tautoty , &
        1056             :            cam_in%shf, cam_in%cflx(:,1), obklen  , ustar       , &
        1057             :            kvm       , kvh      , kvq     , cgh      , cgs     , &
        1058     1489176 :            ri        )
        1059             :
        1060     1489176 :       call pbuf_get_field(pbuf, clubbtop_idx, clubbtop)
        1061             :       !
        1062             :       ! zero out HB where CLUBB is active
        1063             :       !
        1064    24865776 :       do i=1,ncol
        1065   839030941 :         do k=clubbtop(i),pverp
        1066   814165165 :           kvm(i,k) = 0.0_r8
        1067   814165165 :           kvh(i,k) = 0.0_r8
        1068   814165165 :           kvq(i,k) = 0.0_r8
        1069   814165165 :           cgs(i,k) = 0.0_r8
        1070   837541765 :           cgh(i,k) = 0.0_r8
        1071             :         end do
        1072             :       end do
        1073             :
        1074     1489176 :       call outfld( 'HB_ri',          ri,         pcols,   lchnk )
        1075             :     else
        1076             :       ! CLUBB has only a bare-bones placeholder here. If using CLUBB, the
        1077             :       ! PBL diffusion will happen before coupling, so vertical_diffusion
        1078             :       ! is only handling other things, e.g. some boundary conditions, tms,
        1079             :       ! and molecular diffusion.
        1080             :
        1081           0 :       call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol))
        1082             :
        1083           0 :       call calc_ustar( ncol, state%t(:ncol,pver), state%pmid(:ncol,pver), &
        1084           0 :            cam_in%wsx(:ncol), cam_in%wsy(:ncol), rrho(:ncol), ustar(:ncol))
        1085             :       ! Use actual qflux, not lhf/latvap as was done previously
        1086           0 :       call calc_obklen( ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), &
        1087             :            cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol),  &
        1088           0 :            khfs(:ncol), kqfs(:ncol), kbfs(:ncol), obklen(:ncol))
        1089             :       ! These tendencies all applied elsewhere.
        1090           0 :       kvm = 0._r8
        1091           0 :       kvh = 0._r8
        1092           0 :       kvq = 0._r8
        1093             :       ! Not defined since PBL is not actually running here.
        1094           0 :       cgh = 0._r8
        1095           0 :       cgs = 0._r8
        1096             :     end if
        1097             :   end select
        1098             :
        1099     1489176 :   call outfld( 'ustar',   ustar(:), pcols, lchnk )
        1100     1489176 :   call outfld( 'obklen', obklen(:), pcols, lchnk )
        1101             :   !
        1102             :   ! add sponge layer vertical diffusion
        1103             :   !
        1104     1489176 :   if (allocated(kvm_sponge)) then
        1105     7445880 :      do k=1,size(kvm_sponge(:),1)
        1106   100952280 :         kvm(:ncol,1) = kvm(:ncol,1)+kvm_sponge(k)
        1107             :      end do
        1108             :   end if
        1109             :
        1110             :   ! kvh (in pbuf) is used by other physics parameterizations, and as an initial guess in compute_eddy_diff
        1111             :   ! on the next timestep.  It is not updated by the compute_vdiff call below.
        1112     1489176 :   call pbuf_set_field(pbuf, kvh_idx, kvh)
        1113             :
        1114             :   ! kvm (in pbuf) is only used as an initial guess in compute_eddy_diff on the next timestep.
        1115             :   ! The contributions for molecular diffusion made to kvm by the call to compute_vdiff below
        1116             :   ! are not included in the pbuf as these are not needed in the initial guess by compute_eddy_diff.
        1117     1489176 :   call pbuf_set_field(pbuf, kvm_idx, kvm)
        1118             :
        1119             :   !------------------------------------ !
        1120             :   !    Application of diffusivities     !
        1121             :   !------------------------------------ !
        1122             :
        1123             :   ! Set arrays from input state.
        1124 25455558960 :   q_tmp(:ncol,:,:) = state%q(:ncol,:,:)
        1125  2314006344 :   s_tmp(:ncol,:) = state%s(:ncol,:)
        1126  2314006344 :   u_tmp(:ncol,:) = state%u(:ncol,:)
        1127  2314006344 :   v_tmp(:ncol,:) = state%v(:ncol,:)
        1128             :
        1129             :   !------------------------------------------------------ !
        1130             :   ! Write profile output before applying diffusion scheme !
        1131             :   !------------------------------------------------------ !
        1132             :
        1133     1489176 :   if (.not. do_pbl_diags) then
        1134           0 :      sl_prePBL(:ncol,:pver)  = s_tmp(:ncol,:) -   latvap * q_tmp(:ncol,:,ixcldliq) &
        1135           0 :           - ( latvap + latice) * q_tmp(:ncol,:,ixcldice)
        1136             :      qt_prePBL(:ncol,:pver)  = q_tmp(:ncol,:,1) + q_tmp(:ncol,:,ixcldliq) &
        1137           0 :           + q_tmp(:ncol,:,ixcldice)
        1138           0 :      slv_prePBL(:ncol,:pver) = sl_prePBL(:ncol,:pver) * ( 1._r8 + zvir*qt_prePBL(:ncol,:pver) )
        1139             :
        1140           0 :      do k = 1, pver
        1141           0 :         call qsat(state%t(1:ncol,k), state%pmid(1:ncol,k), tem2(1:ncol,k), ftem(1:ncol,k), ncol)
        1142             :      end do
        1143           0 :      ftem_prePBL(:ncol,:) = state%q(:ncol,:,1)/ftem(:ncol,:)*100._r8
        1144             :
        1145           0 :      call outfld( 'qt_pre_PBL   ', qt_prePBL,                 pcols, lchnk )
        1146           0 :      call outfld( 'sl_pre_PBL   ', sl_prePBL,                 pcols, lchnk )
        1147           0 :      call outfld( 'slv_pre_PBL  ', slv_prePBL,                pcols, lchnk )
        1148           0 :      call outfld( 'u_pre_PBL    ', state%u,                   pcols, lchnk )
        1149           0 :      call outfld( 'v_pre_PBL    ', state%v,                   pcols, lchnk )
        1150           0 :      call outfld( 'qv_pre_PBL   ', state%q(:,:,1),            pcols, lchnk )
        1151           0 :      call outfld( 'ql_pre_PBL   ', state%q(:,:,ixcldliq),     pcols, lchnk )
        1152           0 :      call outfld( 'qi_pre_PBL   ', state%q(:,:,ixcldice),     pcols, lchnk )
        1153           0 :      call outfld( 't_pre_PBL    ', state%t,                   pcols, lchnk )
        1154           0 :      call outfld( 'rh_pre_PBL   ', ftem_prePBL,               pcols, lchnk )
        1155             :
        1156             :   end if
        1157             :
        1158             :   ! --------------------------------------------------------------------------------- !
        1159             :   ! Call the diffusivity solver and solve diffusion equation                          !
        1160             :   ! The final two arguments are optional function references to                       !
        1161             :   ! constituent-independent and constituent-dependent moleculuar diffusivity routines !
        1162             :   ! --------------------------------------------------------------------------------- !
        1163             :
        1164             :   ! Modification : We may need to output 'tautotx_im,tautoty_im' from below 'compute_vdiff' and
        1165             :   !                separately print out as diagnostic output, because these are different from
        1166             :   !                the explicit 'tautotx, tautoty' computed above.
        1167             :   ! Note that the output 'tauresx,tauresy' from below subroutines are fully implicit ones.
        1168             :
        1169     1489176 :   call pbuf_get_field(pbuf, kvt_idx, kvt)
        1170             :
        1171     1489176 :   if (do_molec_diff .and. .not. waccmx_mode) then
        1172             :      ! Top boundary condition for dry static energy
        1173           0 :      dse_top(:ncol) = cpairv(:ncol,1,lchnk) * tint(:ncol,1) + &
        1174           0 :           gravit * state%zi(:ncol,1)
        1175             :   else
        1176    24865776 :      dse_top(:ncol) = 0._r8
        1177             :   end if
        1178             :
        1179     1489176 :   select case (eddy_scheme)
        1180             :   case ('CLUBB_SGS')
        1181             :      ! CLUBB applies some fluxes itself, but we still want constituent
        1182             :      ! fluxes applied here (except water vapor).
        1183     1489176 :      taux = 0._r8
        1184     1489176 :      tauy = 0._r8
        1185     1489176 :      shflux = 0._r8
        1186    25315992 :      cflux(:,1) = 0._r8
        1187     1489176 :      if (cam_physpkg_is("cam7")) then
        1188             :        ! surface fluxes applied in clubb emissions module
        1189   254649096 :        cflux(:,2:) = 0._r8
        1190             :      else
        1191           0 :        cflux(:,2:) = cam_in%cflx(:,2:)
        1192             :      end if
        1193             :   case default
        1194           0 :      taux = cam_in%wsx
        1195           0 :      tauy = cam_in%wsy
        1196           0 :      shflux = cam_in%shf
        1197     1489176 :      cflux = cam_in%cflx
        1198             :   end select
        1199             :
        1200     1489176 :   if( any(fieldlist_wet) ) then
        1201             :
        1202     1489176 :      if (do_molec_diff) then
        1203             :         call compute_molec_diff(state%lchnk, pcols, pver, pcnst, ncol, &
        1204             :              kvm, kvt, tint, rhoi, kq_scal, cnst_mw, &
        1205           0 :              mw_fac, nbot_molec)
        1206             :      end if
        1207             :
        1208             :      call compute_vdiff( state%lchnk   ,                                                                     &
        1209             :           pcols         , pver               , pcnst        , ncol          , tint          , &
        1210             :           p    , state%t      , rhoi, ztodt         , taux          , &
        1211             :           tauy          , shflux             , cflux        , &
        1212             :           kvh           , kvm                , kvq          , cgs           , cgh           , &
        1213             :           state%zi      , ksrftms            , dragblj      , &
        1214             :           qmincg       , fieldlist_wet , fieldlist_molec,&
        1215             :           u_tmp         , v_tmp              , q_tmp        , s_tmp         ,                 &
        1216             :           tautmsx       , tautmsy            , dtk          , topflx        , errstring     , &
        1217             :           tauresx       , tauresy            , 1            , cpairv(:,:,state%lchnk), dse_top, &
        1218             :           do_molec_diff, waccmx_mode, &
        1219             :           vd_lu_qdecomp, &
        1220             :           ubc_mmr, ubc_flux, kvt, state%pmid, &
        1221             :           cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx, nbot_molec, &
        1222     1489176 :           kq_scal, mw_fac)
        1223             :
        1224             :      call handle_errmsg(errstring, subname="compute_vdiff", &
        1225     1489176 :           extra_msg="Error in fieldlist_wet call from vertical_diffusion.")
        1226             :
        1227             :   end if
        1228             :
        1229     1489176 :   if( any( fieldlist_dry ) ) then
        1230             :
        1231           0 :      if( do_molec_diff ) then
        1232             :         ! kvm is unused in the output here (since it was assigned
        1233             :         ! above), so we use a temp kvm for the inout argument, and
        1234             :         ! ignore the value output by compute_molec_diff.
        1235           0 :         kvm_temp = kvm
        1236             :         call compute_molec_diff(state%lchnk, pcols, pver, pcnst, ncol, &
        1237             :              kvm_temp, kvt, tint, rhoi_dry, kq_scal, cnst_mw, &
        1238           0 :              mw_fac, nbot_molec)
        1239             :      end if
        1240             :
        1241             :      call compute_vdiff( state%lchnk   ,                                                                     &
        1242             :           pcols         , pver               , pcnst        , ncol          , tint          , &
        1243             :           p_dry , state%t      , rhoi_dry,  ztodt         , taux          , &
        1244             :           tauy          , shflux             , cflux        , &
        1245             :           kvh           , kvm                , kvq          , cgs           , cgh           , &
        1246             :           state%zi      , ksrftms            , dragblj      , &
        1247             :           qmincg       , fieldlist_dry , fieldlist_molec,&
        1248             :           u_tmp         , v_tmp              , q_tmp        , s_tmp         ,                 &
        1249             :           tautmsx_temp  , tautmsy_temp       , dtk_temp     , topflx_temp   , errstring     , &
        1250             :           tauresx       , tauresy            , 1            , cpairv(:,:,state%lchnk), dse_top, &
        1251             :           do_molec_diff , waccmx_mode, &
        1252             :           vd_lu_qdecomp, &
        1253             :           ubc_mmr, ubc_flux, kvt, state%pmiddry, &
        1254             :           cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx, nbot_molec, &
        1255           0 :           kq_scal, mw_fac)
        1256             :
        1257             :      call handle_errmsg(errstring, subname="compute_vdiff", &
        1258           0 :           extra_msg="Error in fieldlist_dry call from vertical_diffusion.")
        1259             :
        1260             :   end if
        1261             :
        1262     1489176 :   if (prog_modal_aero) then
        1263             :
        1264             :      ! Modal aerosol species not diffused, so just add the explicit surface fluxes to the
        1265             :      ! lowest layer.  **NOTE** This code assumes wet mmr.
        1266             :
        1267           0 :      tmp1(:ncol) = ztodt * gravit * state%rpdel(:ncol,pver)
        1268           0 :      do m = 1, pmam_ncnst
        1269           0 :         l = pmam_cnst_idx(m)
        1270           0 :         q_tmp(:ncol,pver,l) = q_tmp(:ncol,pver,l) + tmp1(:ncol) * cflux(:ncol,l)
        1271             :      enddo
        1272             :   end if
        1273             :
        1274             :   ! -------------------------------------------------------- !
        1275             :   ! Diagnostics and output writing after applying PBL scheme !
        1276             :   ! -------------------------------------------------------- !
        1277             :
        1278     1489176 :   if (.not. do_pbl_diags) then
        1279             :
        1280           0 :      sl(:ncol,:pver)  = s_tmp(:ncol,:) -   latvap           * q_tmp(:ncol,:,ixcldliq) &
        1281           0 :           - ( latvap + latice) * q_tmp(:ncol,:,ixcldice)
        1282             :      qt(:ncol,:pver)  = q_tmp(:ncol,:,1) + q_tmp(:ncol,:,ixcldliq) &
        1283           0 :           + q_tmp(:ncol,:,ixcldice)
        1284           0 :      slv(:ncol,:pver) = sl(:ncol,:pver) * ( 1._r8 + zvir*qt(:ncol,:pver) )
        1285             :
        1286           0 :      slflx(:ncol,1) = 0._r8
        1287           0 :      qtflx(:ncol,1) = 0._r8
        1288           0 :      uflx(:ncol,1)  = 0._r8
        1289           0 :      vflx(:ncol,1)  = 0._r8
        1290             :
        1291           0 :      slflx_cg(:ncol,1) = 0._r8
        1292           0 :      qtflx_cg(:ncol,1) = 0._r8
        1293           0 :      uflx_cg(:ncol,1)  = 0._r8
        1294           0 :      vflx_cg(:ncol,1)  = 0._r8
        1295             :
        1296           0 :      do k = 2, pver
        1297           0 :         do i = 1, ncol
        1298           0 :            rhoair     = state%pint(i,k) / ( rair * ( ( 0.5_r8*(slv(i,k)+slv(i,k-1)) - gravit*state%zi(i,k))/cpair ) )
        1299           0 :            slflx(i,k) = kvh(i,k) * &
        1300           0 :                 ( - rhoair*(sl(i,k-1)-sl(i,k))/(state%zm(i,k-1)-state%zm(i,k)) &
        1301           0 :                 + cgh(i,k) )
        1302             :            qtflx(i,k) = kvh(i,k) * &
        1303           0 :                 ( - rhoair*(qt(i,k-1)-qt(i,k))/(state%zm(i,k-1)-state%zm(i,k)) &
        1304           0 :                 + rhoair*(cam_in%cflx(i,1)+cam_in%cflx(i,ixcldliq)+cam_in%cflx(i,ixcldice))*cgs(i,k) )
        1305             :            uflx(i,k)  = kvm(i,k) * &
        1306           0 :                 ( - rhoair*(u_tmp(i,k-1)-u_tmp(i,k))/(state%zm(i,k-1)-state%zm(i,k)))
        1307             :            vflx(i,k)  = kvm(i,k) * &
        1308           0 :                 ( - rhoair*(v_tmp(i,k-1)-v_tmp(i,k))/(state%zm(i,k-1)-state%zm(i,k)))
        1309           0 :            slflx_cg(i,k) = kvh(i,k) * cgh(i,k)
        1310             :            qtflx_cg(i,k) = kvh(i,k) * rhoair * ( cam_in%cflx(i,1) + &
        1311           0 :                 cam_in%cflx(i,ixcldliq) + cam_in%cflx(i,ixcldice) ) * cgs(i,k)
        1312           0 :            uflx_cg(i,k)  = 0._r8
        1313           0 :            vflx_cg(i,k)  = 0._r8
        1314             :         end do
        1315             :      end do
        1316             :
        1317             :      ! Modification : I should check whether slflx(:ncol,pverp) is correctly computed.
        1318             :      !                Note also that 'tautotx' is explicit total stress, different from
        1319             :      !                the ones that have been actually added into the atmosphere.
        1320             :
        1321           0 :      slflx(:ncol,pverp) = cam_in%shf(:ncol)
        1322           0 :      qtflx(:ncol,pverp) = cam_in%cflx(:ncol,1)
        1323           0 :      uflx(:ncol,pverp)  = tautotx(:ncol)
        1324           0 :      vflx(:ncol,pverp)  = tautoty(:ncol)
        1325             :
        1326           0 :      slflx_cg(:ncol,pverp) = 0._r8
        1327           0 :      qtflx_cg(:ncol,pverp) = 0._r8
        1328           0 :      uflx_cg(:ncol,pverp)  = 0._r8
        1329           0 :      vflx_cg(:ncol,pverp)  = 0._r8
        1330             :
        1331           0 :      if (trim(shallow_scheme) == 'UNICON') then
        1332           0 :         call pbuf_get_field(pbuf, qtl_flx_idx,  qtl_flx)
        1333           0 :         call pbuf_get_field(pbuf, qti_flx_idx,  qti_flx)
        1334           0 :         qtl_flx(:ncol,1) = 0._r8
        1335           0 :         qti_flx(:ncol,1) = 0._r8
        1336           0 :         do k = 2, pver
        1337           0 :            do i = 1, ncol
        1338             :               ! For use in the cloud macrophysics
        1339             :               ! Note that density is not added here. Also, only consider local transport term.
        1340           0 :               qtl_flx(i,k) = - kvh(i,k)*(q_tmp(i,k-1,1)-q_tmp(i,k,1)+q_tmp(i,k-1,ixcldliq)-q_tmp(i,k,ixcldliq))/&
        1341           0 :                    (state%zm(i,k-1)-state%zm(i,k))
        1342           0 :               qti_flx(i,k) = - kvh(i,k)*(q_tmp(i,k-1,1)-q_tmp(i,k,1)+q_tmp(i,k-1,ixcldice)-q_tmp(i,k,ixcldice))/&
        1343           0 :                    (state%zm(i,k-1)-state%zm(i,k))
        1344             :            end do
        1345             :         end do
        1346           0 :         do i = 1, ncol
        1347           0 :            rhoair = state%pint(i,pverp)/(rair*((slv(i,pver)-gravit*state%zi(i,pverp))/cpair))
        1348           0 :            qtl_flx(i,pverp) = cam_in%cflx(i,1)/rhoair
        1349           0 :            qti_flx(i,pverp) = cam_in%cflx(i,1)/rhoair
        1350             :         end do
        1351             :      end if
        1352             :
        1353             :   end if
        1354             :
        1355             :   ! --------------------------------------------------------------- !
        1356             :   ! Convert the new profiles into vertical diffusion tendencies.    !
        1357             :   ! Convert KE dissipative heat change into "temperature" tendency. !
        1358             :   ! --------------------------------------------------------------- !
        1359             :   ! All variables are modified by vertical diffusion
        1360             :
        1361    17870112 :   lq(:) = .TRUE.
        1362             :   call physics_ptend_init(ptend,state%psetcols, "vertical diffusion", &
        1363     1489176 :        ls=.true., lu=.true., lv=.true., lq=lq)
        1364             :
        1365  2314006344 :   ptend%s(:ncol,:)       = ( s_tmp(:ncol,:) - state%s(:ncol,:) ) * rztodt
        1366  2314006344 :   ptend%u(:ncol,:)       = ( u_tmp(:ncol,:) - state%u(:ncol,:) ) * rztodt
        1367  2314006344 :   ptend%v(:ncol,:)       = ( v_tmp(:ncol,:) - state%v(:ncol,:) ) * rztodt
        1368 25455558960 :   ptend%q(:ncol,:pver,:) = ( q_tmp(:ncol,:pver,:) - state%q(:ncol,:pver,:) ) * rztodt
        1369             :
        1370             :   ! Convert tendencies of dry constituents to dry basis.
        1371    17870112 :   do m = 1,pcnst
        1372    17870112 :      if (cnst_type(m).eq.'dry') then
        1373           0 :         ptend%q(:ncol,:pver,m) = ptend%q(:ncol,:pver,m)*state%pdel(:ncol,:pver)/state%pdeldry(:ncol,:pver)
        1374             :      endif
        1375             :   end do
        1376             :   ! convert wet mmr back to dry before conservation check
        1377     1489176 :   call set_wet_to_dry(state, convert_cnst_type='dry')
        1378             :
        1379     1489176 :   if (.not. do_pbl_diags) then
        1380           0 :      slten(:ncol,:)         = ( sl(:ncol,:) - sl_prePBL(:ncol,:) ) * rztodt
        1381           0 :      qtten(:ncol,:)         = ( qt(:ncol,:) - qt_prePBL(:ncol,:) ) * rztodt
        1382             :   end if
        1383             :
        1384             :   ! ------------------------------------------------------------ !
        1385             :   ! In order to perform 'pseudo-conservative variable diffusion' !
        1386             :   ! perform the following two stages:                            !
        1387             :   !                                                              !
        1388             :   ! I.  Re-set (1) 'qvten' by 'qtten', and 'qlten = qiten = 0'   !
        1389             :   !            (2) 'sten'  by 'slten', and                       !
        1390             :   !            (3) 'qlten = qiten = 0'                           !
        1391             :   !                                                              !
        1392             :   ! II. Apply 'positive_moisture'                                !
        1393             :   !                                                              !
        1394             :   ! ------------------------------------------------------------ !
        1395             :
        1396     1489176 :     if( (eddy_scheme .eq. 'diag_TKE' .or. eddy_scheme .eq. 'SPCAM_m2005') .and. do_pseudocon_diff ) then
        1397             :
        1398           0 :      ptend%q(:ncol,:pver,1) = qtten(:ncol,:pver)
        1399           0 :      ptend%s(:ncol,:pver)   = slten(:ncol,:pver)
        1400           0 :      ptend%q(:ncol,:pver,ixcldliq) = 0._r8
        1401           0 :      ptend%q(:ncol,:pver,ixcldice) = 0._r8
        1402           0 :      if (ixnumliq > 0) ptend%q(:ncol,:pver,ixnumliq) = 0._r8
        1403           0 :      if (ixnumice > 0) ptend%q(:ncol,:pver,ixnumice) = 0._r8
        1404             :
        1405           0 :      do i = 1, ncol
        1406           0 :         do k = 1, pver
        1407           0 :            qv_pro(i,k) = state%q(i,k,1)        + ptend%q(i,k,1)             * ztodt
        1408           0 :            ql_pro(i,k) = state%q(i,k,ixcldliq) + ptend%q(i,k,ixcldliq)      * ztodt
        1409           0 :            qi_pro(i,k) = state%q(i,k,ixcldice) + ptend%q(i,k,ixcldice)      * ztodt
        1410           0 :            s_pro(i,k)  = state%s(i,k)          + ptend%s(i,k)               * ztodt
        1411           0 :            t_pro(i,k)  = state%t(i,k)          + (1._r8/cpair)*ptend%s(i,k) * ztodt
        1412             :         end do
        1413             :      end do
        1414             :
        1415           0 :      call positive_moisture( cpair, latvap, latvap+latice, ncol, pver, ztodt, qmin(1), qmin(ixcldliq), qmin(ixcldice),    &
        1416           0 :           state%pdel(:ncol,pver:1:-1), qv_pro(:ncol,pver:1:-1), ql_pro(:ncol,pver:1:-1), &
        1417           0 :           qi_pro(:ncol,pver:1:-1), t_pro(:ncol,pver:1:-1), s_pro(:ncol,pver:1:-1),       &
        1418           0 :           ptend%q(:ncol,pver:1:-1,1), ptend%q(:ncol,pver:1:-1,ixcldliq),                 &
        1419           0 :           ptend%q(:ncol,pver:1:-1,ixcldice), ptend%s(:ncol,pver:1:-1) )
        1420             :
        1421             :   end if
        1422             :
        1423             :   ! ----------------------------------------------------------------- !
        1424             :   ! Re-calculate diagnostic output variables after vertical diffusion !
        1425             :   ! ----------------------------------------------------------------- !
        1426             :
        1427     1489176 :   if (.not. do_pbl_diags) then
        1428             :
        1429           0 :      qv_aft_PBL(:ncol,:pver)  =   state%q(:ncol,:pver,1)         + ptend%q(:ncol,:pver,1)        * ztodt
        1430           0 :      ql_aft_PBL(:ncol,:pver)  =   state%q(:ncol,:pver,ixcldliq)  + ptend%q(:ncol,:pver,ixcldliq) * ztodt
        1431           0 :      qi_aft_PBL(:ncol,:pver)  =   state%q(:ncol,:pver,ixcldice)  + ptend%q(:ncol,:pver,ixcldice) * ztodt
        1432           0 :      s_aft_PBL(:ncol,:pver)   =   state%s(:ncol,:pver)           + ptend%s(:ncol,:pver)          * ztodt
        1433           0 :      t_aftPBL(:ncol,:pver)    = ( s_aft_PBL(:ncol,:pver) - gravit*state%zm(:ncol,:pver) ) / cpair
        1434             :
        1435           0 :      u_aft_PBL(:ncol,:pver)   =  state%u(:ncol,:pver)          + ptend%u(:ncol,:pver)            * ztodt
        1436           0 :      v_aft_PBL(:ncol,:pver)   =  state%v(:ncol,:pver)          + ptend%v(:ncol,:pver)            * ztodt
        1437             :
        1438           0 :      do k = 1, pver
        1439           0 :         call qsat(t_aftPBL(1:ncol,k), state%pmid(1:ncol,k), tem2(1:ncol,k), ftem(1:ncol,k), ncol)
        1440             :      end do
        1441           0 :      ftem_aftPBL(:ncol,:pver) = qv_aft_PBL(:ncol,:pver) / ftem(:ncol,:pver) * 100._r8
        1442             :
        1443           0 :      tten(:ncol,:pver)        = ( t_aftPBL(:ncol,:pver)    - state%t(:ncol,:pver) )              * rztodt
        1444           0 :      rhten(:ncol,:pver)       = ( ftem_aftPBL(:ncol,:pver) - ftem_prePBL(:ncol,:pver) )          * rztodt
        1445             :
        1446             :   end if
        1447             :
        1448             :
        1449             :   ! -------------------------------------------------------------- !
        1450             :   ! mass conservation check.........
        1451             :   ! -------------------------------------------------------------- !
        1452     1489176 :   if (diff_cnsrv_mass_check) then
        1453             :
        1454             :      ! Conservation check
        1455           0 :      do m = 1, pcnst
        1456           0 :         fixed_ubc: if ((.not.cnst_fixed_ubc(m)).and.(.not.cnst_fixed_ubflx(m))) then
        1457           0 :            col_loop: do i = 1, ncol
        1458           0 :               sum1 = 0._r8
        1459           0 :               sum2 = 0._r8
        1460           0 :               sum3 = 0._r8
        1461           0 :               do k = 1, pver
        1462           0 :                  if(cnst_get_type_byind(m).eq.'wet') then
        1463           0 :                     pdelx = state%pdel(i,k)
        1464             :                  else
        1465           0 :                     pdelx = state%pdeldry(i,k)
        1466             :                  endif
        1467           0 :                  sum1 = sum1 + state%q(i,k,m)*pdelx/gravit                          ! total column
        1468           0 :                  sum2 = sum2 +(state%q(i,k,m)+ptend%q(i,k,m)*ztodt)*pdelx/ gravit   ! total column after tendancy is applied
        1469           0 :                  sum3 = sum3 +(               ptend%q(i,k,m)*ztodt)*pdelx/ gravit   ! rate of change in column
        1470             :               enddo
        1471           0 :               sum1 = sum1 + (cam_in%cflx(i,m) * ztodt) ! add in surface flux (kg/m2)
        1472           0 :               sflx = (cam_in%cflx(i,m) * ztodt)
        1473           0 :               if (sum1>1.e-36_r8) then
        1474           0 :                  if( abs((sum2-sum1)/sum1) .gt. 1.e-12_r8  ) then
        1475           0 :                     nstep = get_nstep()
        1476             :                     write(iulog,'(a,a8,a,I4,2f8.3,5e25.16)') &
        1477           0 :                          'MASSCHECK vert diff : nstep,lon,lat,mass1,mass2,sum3,sflx,rel-diff : ', &
        1478           0 :                          trim(cnst_name(m)), ' : ', nstep, state%lon(i)*180._r8/pi, state%lat(i)*180._r8/pi, &
        1479           0 :                          sum1, sum2, sum3, sflx, abs(sum2-sum1)/sum1
        1480             : !xxx                    call endrun('vertical_diffusion_tend : mass not conserved' )
        1481             :                  endif
        1482             :               endif
        1483             :            enddo col_loop
        1484             :         endif fixed_ubc
        1485             :      enddo
        1486             :   endif
        1487             :
        1488             :   ! -------------------------------------------------------------- !
        1489             :   ! Writing state variables after PBL scheme for detailed analysis !
        1490             :   ! -------------------------------------------------------------- !
        1491             :
        1492     1489176 :   if (.not. do_pbl_diags) then
        1493             :
        1494           0 :      call outfld( 'sl_aft_PBL'   , sl,                        pcols, lchnk )
        1495           0 :      call outfld( 'qt_aft_PBL'   , qt,                        pcols, lchnk )
        1496           0 :      call outfld( 'slv_aft_PBL'  , slv,                       pcols, lchnk )
        1497           0 :      call outfld( 'u_aft_PBL'    , u_aft_PBL,                 pcols, lchnk )
        1498           0 :      call outfld( 'v_aft_PBL'    , v_aft_PBL,                 pcols, lchnk )
        1499           0 :      call outfld( 'qv_aft_PBL'   , qv_aft_PBL,                pcols, lchnk )
        1500           0 :      call outfld( 'ql_aft_PBL'   , ql_aft_PBL,                pcols, lchnk )
        1501           0 :      call outfld( 'qi_aft_PBL'   , qi_aft_PBL,                pcols, lchnk )
        1502           0 :      call outfld( 't_aft_PBL '   , t_aftPBL,                  pcols, lchnk )
        1503           0 :      call outfld( 'rh_aft_PBL'   , ftem_aftPBL,               pcols, lchnk )
        1504           0 :      call outfld( 'slflx_PBL'    , slflx,                     pcols, lchnk )
        1505           0 :      call outfld( 'qtflx_PBL'    , qtflx,                     pcols, lchnk )
        1506           0 :      call outfld( 'uflx_PBL'     , uflx,                      pcols, lchnk )
        1507           0 :      call outfld( 'vflx_PBL'     , vflx,                      pcols, lchnk )
        1508           0 :      call outfld( 'slflx_cg_PBL' , slflx_cg,                  pcols, lchnk )
        1509           0 :      call outfld( 'qtflx_cg_PBL' , qtflx_cg,                  pcols, lchnk )
        1510           0 :      call outfld( 'uflx_cg_PBL'  , uflx_cg,                   pcols, lchnk )
        1511           0 :      call outfld( 'vflx_cg_PBL'  , vflx_cg,                   pcols, lchnk )
        1512           0 :      call outfld( 'slten_PBL'    , slten,                     pcols, lchnk )
        1513           0 :      call outfld( 'qtten_PBL'    , qtten,                     pcols, lchnk )
        1514           0 :      call outfld( 'uten_PBL'     , ptend%u,                   pcols, lchnk )
        1515           0 :      call outfld( 'vten_PBL'     , ptend%v,                   pcols, lchnk )
        1516           0 :      call outfld( 'qvten_PBL'    , ptend%q(:,:,1),            pcols, lchnk )
        1517           0 :      call outfld( 'qlten_PBL'    , ptend%q(:,:,ixcldliq),     pcols, lchnk )
        1518           0 :      call outfld( 'qiten_PBL'    , ptend%q(:,:,ixcldice),     pcols, lchnk )
        1519           0 :      call outfld( 'tten_PBL'     , tten,                      pcols, lchnk )
        1520           0 :      call outfld( 'rhten_PBL'    , rhten,                     pcols, lchnk )
        1521             :
        1522             :   end if
        1523             :
        1524             :   ! ------------------------------------------- !
        1525             :   ! Writing the other standard output variables !
        1526             :   ! ------------------------------------------- !
        1527             :
        1528     1489176 :   if (.not. do_pbl_diags) then
        1529           0 :      call outfld( 'QT'           , qt,                        pcols, lchnk )
        1530           0 :      call outfld( 'SL'           , sl,                        pcols, lchnk )
        1531           0 :      call outfld( 'SLV'          , slv,                       pcols, lchnk )
        1532           0 :      call outfld( 'SLFLX'        , slflx,                     pcols, lchnk )
        1533           0 :      call outfld( 'QTFLX'        , qtflx,                     pcols, lchnk )
        1534           0 :      call outfld( 'UFLX'         , uflx,                      pcols, lchnk )
        1535           0 :      call outfld( 'VFLX'         , vflx,                      pcols, lchnk )
        1536           0 :      call outfld( 'TKE'          , tke,                       pcols, lchnk )
        1537             :
        1538           0 :      call outfld( 'PBLH'         , pblh,                      pcols, lchnk )
        1539           0 :      call outfld( 'TPERT'        , tpert,                     pcols, lchnk )
        1540           0 :      call outfld( 'QPERT'        , qpert,                     pcols, lchnk )
        1541             :   end if
        1542     1489176 :   call outfld( 'USTAR'        , ustar,                     pcols, lchnk )
        1543     1489176 :   call outfld( 'KVH'          , kvh,                       pcols, lchnk )
        1544     1489176 :   call outfld( 'KVT'          , kvt,                       pcols, lchnk )
        1545     1489176 :   call outfld( 'KVM'          , kvm,                       pcols, lchnk )
        1546     1489176 :   call outfld( 'CGS'          , cgs,                       pcols, lchnk )
        1547  2314006344 :   dtk(:ncol,:) = dtk(:ncol,:) / cpair / ztodt      ! Normalize heating for history
        1548     1489176 :   call outfld( 'DTVKE'        , dtk,                       pcols, lchnk )
        1549  2314006344 :   dtk(:ncol,:) = ptend%s(:ncol,:) / cpair          ! Normalize heating for history using dtk
        1550     1489176 :   call outfld( 'DTV'          , dtk,                       pcols, lchnk )
        1551     1489176 :   call outfld( 'DUV'          , ptend%u,                   pcols, lchnk )
        1552     1489176 :   call outfld( 'DVV'          , ptend%v,                   pcols, lchnk )
        1553    17870112 :   do m = 1, pcnst
        1554    17870112 :      call outfld( vdiffnam(m) , ptend%q(1,1,m),            pcols, lchnk )
        1555             :   end do
        1556     1489176 :   if( do_molec_diff ) then
        1557           0 :      call outfld( 'TTPXMLC'  , topflx,                    pcols, lchnk )
        1558             :   end if
        1559             :
        1560     1489176 :   call p%finalize()
        1561     1489176 :   call p_dry%finalize()
        1562             :
        1563     2978352 : end subroutine vertical_diffusion_tend
        1564             :
        1565             : ! =============================================================================== !
        1566             : !                                                                                 !
        1567             : ! =============================================================================== !
        1568             :
        1569           0 : subroutine positive_moisture( cp, xlv, xls, ncol, mkx, dt, qvmin, qlmin, qimin, &
        1570           0 :      dp, qv, ql, qi, t, s, qvten, qlten, qiten, sten )
        1571             :   ! ------------------------------------------------------------------------------- !
        1572             :   ! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer,         !
        1573             :   ! force them to be larger than minimum value by (1) condensating water vapor      !
        1574             :   ! into liquid or ice, and (2) by transporting water vapor from the very lower     !
        1575             :   ! layer. '2._r8' is multiplied to the minimum values for safety.                  !
        1576             :   ! Update final state variables and tendencies associated with this correction.    !
        1577             :   ! If any condensation happens, update (s,t) too.                                  !
        1578             :   ! Note that (qv,ql,qi,t,s) are final state variables after applying corresponding !
        1579             :   ! input tendencies.                                                               !
        1580             :   ! Be careful the order of k : '1': near-surface layer, 'mkx' : top layer          !
        1581             :   ! ------------------------------------------------------------------------------- !
        1582             :   implicit none
        1583             :   integer,  intent(in)     :: ncol, mkx
        1584             :   real(r8), intent(in)     :: cp, xlv, xls
        1585             :   real(r8), intent(in)     :: dt, qvmin, qlmin, qimin
        1586             :   real(r8), intent(in)     :: dp(ncol,mkx)
        1587             :   real(r8), intent(inout)  :: qv(ncol,mkx), ql(ncol,mkx), qi(ncol,mkx), t(ncol,mkx), s(ncol,mkx)
        1588             :   real(r8), intent(inout)  :: qvten(ncol,mkx), qlten(ncol,mkx), qiten(ncol,mkx), sten(ncol,mkx)
        1589             :   integer   i, k
        1590             :   real(r8)  dql, dqi, dqv, sum, aa, dum
        1591             :
        1592             :   ! Modification : I should check whether this is exactly same as the one used in
        1593             :   !                shallow convection and cloud macrophysics.
        1594             :
        1595           0 :   do i = 1, ncol
        1596           0 :      do k = mkx, 1, -1    ! From the top to the 1st (lowest) layer from the surface
        1597           0 :         dql        = max(0._r8,1._r8*qlmin-ql(i,k))
        1598           0 :         dqi        = max(0._r8,1._r8*qimin-qi(i,k))
        1599           0 :         qlten(i,k) = qlten(i,k) +  dql/dt
        1600           0 :         qiten(i,k) = qiten(i,k) +  dqi/dt
        1601           0 :         qvten(i,k) = qvten(i,k) - (dql+dqi)/dt
        1602           0 :         sten(i,k)  = sten(i,k)  + xlv * (dql/dt) + xls * (dqi/dt)
        1603           0 :         ql(i,k)    = ql(i,k) +  dql
        1604           0 :         qi(i,k)    = qi(i,k) +  dqi
        1605           0 :         qv(i,k)    = qv(i,k) -  dql - dqi
        1606           0 :         s(i,k)     = s(i,k)  +  xlv * dql + xls * dqi
        1607           0 :         t(i,k)     = t(i,k)  + (xlv * dql + xls * dqi)/cp
        1608           0 :         dqv        = max(0._r8,1._r8*qvmin-qv(i,k))
        1609           0 :         qvten(i,k) = qvten(i,k) + dqv/dt
        1610           0 :         qv(i,k)    = qv(i,k)    + dqv
        1611           0 :         if( k .ne. 1 ) then
        1612           0 :            qv(i,k-1)    = qv(i,k-1)    - dqv*dp(i,k)/dp(i,k-1)
        1613           0 :            qvten(i,k-1) = qvten(i,k-1) - dqv*dp(i,k)/dp(i,k-1)/dt
        1614             :         endif
        1615           0 :         qv(i,k) = max(qv(i,k),qvmin)
        1616           0 :         ql(i,k) = max(ql(i,k),qlmin)
        1617           0 :         qi(i,k) = max(qi(i,k),qimin)
        1618             :      end do
        1619             :      ! Extra moisture used to satisfy 'qv(i,1)=qvmin' is proportionally
        1620             :      ! extracted from all the layers that has 'qv > 2*qvmin'. This fully
        1621             :      ! preserves column moisture.
        1622           0 :      if( dqv .gt. 1.e-20_r8 ) then
        1623             :         sum = 0._r8
        1624           0 :         do k = 1, mkx
        1625           0 :            if( qv(i,k) .gt. 2._r8*qvmin ) sum = sum + qv(i,k)*dp(i,k)
        1626             :         enddo
        1627           0 :         aa = dqv*dp(i,1)/max(1.e-20_r8,sum)
        1628           0 :         if( aa .lt. 0.5_r8 ) then
        1629           0 :            do k = 1, mkx
        1630           0 :               if( qv(i,k) .gt. 2._r8*qvmin ) then
        1631           0 :                  dum        = aa*qv(i,k)
        1632           0 :                  qv(i,k)    = qv(i,k) - dum
        1633           0 :                  qvten(i,k) = qvten(i,k) - dum/dt
        1634             :               endif
        1635             :            enddo
        1636             :         else
        1637           0 :            write(iulog,*) 'Full positive_moisture is impossible in vertical_diffusion'
        1638             :         endif
        1639             :      endif
        1640             :   end do
        1641           0 :   return
        1642             :
        1643     1489176 : end subroutine positive_moisture
        1644             :
        1645             : end module vertical_diffusion

Generated by: LCOV version 1.14