LCOV - code coverage report
Current view: | top level - physics/cam - diffusion_solver.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 212 | 263 | 80.6 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 6 | 9 | 66.7 %

              Line data    Source code
           1             :
           2             :   module diffusion_solver
           3             :
           4             :   !------------------------------------------------------------------------------------ !
           5             :   ! Module to solve vertical diffusion equations using a tri-diagonal solver.           !
           6             :   ! The module will also apply countergradient fluxes, and apply molecular              !
           7             :   ! diffusion for constituents.                                                         !
           8             :   !                                                                                     !
           9             :   ! Public interfaces :                                                                 !
          10             :   !    init_vdiff       initializes time independent coefficients                       !
          11             :   !    compute_vdiff    solves diffusion equations                                      !
          12             :   !    vdiff_selector   type for storing fields selected to be diffused                 !
          13             :   !    vdiff_select     selects fields to be diffused                                   !
          14             :   !    operator(.not.)  extends .not. to operate on type vdiff_selector                 !
          15             :   !    any              provides functionality of intrinsic any for type vdiff_selector !
          16             :   !                                                                                     !
          17             :   !------------------------------------ Code History ---------------------------------- !
          18             :   ! Initial subroutines :  B. Boville and others, 1991-2004                             !
          19             :   ! Modularization      :  J. McCaa, September 2004                                     !
          20             :   ! Most Recent Code    :  Sungsu Park, Aug. 2006, Dec. 2008, Jan. 2010.                !
          21             :   !------------------------------------------------------------------------------------ !
          22             :
          23             :   implicit none
          24             :   private
          25             :   save
          26             :
          27             :   integer, parameter :: r8 = selected_real_kind(12)      ! 8 byte real
          28             :
          29             :   ! ----------------- !
          30             :   ! Public interfaces !
          31             :   ! ----------------- !
          32             :
          33             :   public init_vdiff                                      ! Initialization
          34             :   public new_fieldlist_vdiff                             ! Returns an empty fieldlist
          35             :   public compute_vdiff                                   ! Full routine
          36             :   public vdiff_selector                                  ! Type for storing fields selected to be diffused
          37             :   public vdiff_select                                    ! Selects fields to be diffused
          38             :   public operator(.not.)                                 ! Extends .not. to operate on type vdiff_selector
          39             :   public any                                             ! Provides functionality of intrinsic any for type vdiff_selector
          40             :
          41             :   ! Below stores logical array of fields to be diffused
          42             :
          43             :   type vdiff_selector
          44             :        private
          45             :        logical, allocatable, dimension(:) :: fields
          46             :   end type vdiff_selector
          47             :
          48             :   ! Below extends .not. to operate on type vdiff_selector
          49             :
          50             :   interface operator(.not.)
          51             :        module procedure not
          52             :   end interface
          53             :
          54             :   ! Below provides functionality of intrinsic any for type vdiff_selector
          55             :
          56             :   interface any
          57             :        module procedure my_any
          58             :   end interface
          59             :
          60             :   ! ------------ !
          61             :   ! Private data !
          62             :   ! ------------ !
          63             :
          64             :   ! Unit number for log output
          65             :   integer :: iulog = -1
          66             :
          67             :   real(r8), private   :: cpair                           ! Specific heat of dry air
          68             :   real(r8), private   :: gravit                          ! Acceleration due to gravity
          69             :   real(r8), private   :: rair                            ! Gas constant for dry air
          70             :
          71             :   logical,  private   :: do_iss                          ! Use implicit turbulent surface stress computation
          72             :
          73             :   ! Parameters used for Turbulent Mountain Stress
          74             :
          75             :   real(r8), parameter :: z0fac   = 0.025_r8              ! Factor determining z_0 from orographic standard deviation
          76             :   real(r8), parameter :: z0max   = 100._r8               ! Max value of z_0 for orography
          77             :   real(r8), parameter :: horomin = 10._r8                ! Min value of subgrid orographic height for mountain stress
          78             :   real(r8), parameter :: dv2min  = 0.01_r8               ! Minimum shear squared
          79             :
          80             :   logical :: am_correction ! logical switch for AM correction
          81             :
          82             :   contains
          83             :
          84             :   ! =============================================================================== !
          85             :   !                                                                                 !
          86             :   ! =============================================================================== !
          87             :
          88        1536 :   subroutine init_vdiff( kind, iulog_in, rair_in, cpair_in, gravit_in, do_iss_in, &
          89        1536 :                          am_correction_in, errstring )
          90             :
          91             :     integer,              intent(in)  :: kind            ! Kind used for reals
          92             :     integer,              intent(in)  :: iulog_in        ! Unit number for log output.
          93             :     real(r8),             intent(in)  :: rair_in         ! Input gas constant for dry air
          94             :     real(r8),             intent(in)  :: cpair_in        ! Input heat capacity for dry air
          95             :     real(r8),             intent(in)  :: gravit_in       ! Input gravitational acceleration
          96             :     logical,              intent(in)  :: do_iss_in       ! Input ISS flag
          97             :     logical,              intent(in)  :: am_correction_in! for angular momentum conservation
          98             :     character(128),       intent(out) :: errstring       ! Output status
          99             :
         100        1536 :     errstring = ''
         101        1536 :     iulog = iulog_in
         102        1536 :     if( kind .ne. r8 ) then
         103           0 :         write(iulog,*) 'KIND of reals passed to init_vdiff -- exiting.'
         104           0 :         errstring = 'init_vdiff'
         105           0 :         return
         106             :     endif
         107             :
         108        1536 :     rair   = rair_in
         109        1536 :     cpair  = cpair_in
         110        1536 :     gravit = gravit_in
         111        1536 :     do_iss = do_iss_in
         112        1536 :     am_correction = am_correction_in
         113             :
         114             :   end subroutine init_vdiff
         115             :
         116             :   ! =============================================================================== !
         117             :   !                                                                                 !
         118             :   ! =============================================================================== !
         119             :
         120        4608 :   type(vdiff_selector) pure function new_fieldlist_vdiff(ncnst)
         121             :
         122             :     integer,              intent(in)  :: ncnst           ! Number of constituents
         123             :
         124       13824 :     allocate( new_fieldlist_vdiff%fields( 3 + ncnst ) )
         125       69120 :     new_fieldlist_vdiff%fields = .false.
         126             :
         127        4608 :   end function new_fieldlist_vdiff
         128             :
         129             :   ! =============================================================================== !
         130             :   !                                                                                 !
         131             :   ! =============================================================================== !
         132             :
         133     1489176 :   subroutine compute_vdiff( lchnk           ,                                                                   &
         134     1489176 :                             pcols           , pver               , ncnst         , ncol         , tint        , &
         135     2978352 :                             p               , t                  , rhoi          , ztodt        , taux        , &
         136     1489176 :                             tauy            , shflx              , cflx          , &
         137     1489176 :                             kvh             , kvm                , kvq           , cgs          , cgh         , &
         138     1489176 :                             zi              , ksrftms            , dragblj       , &
         139     1489176 :                             qmincg          , fieldlist          , fieldlistm    , &
         140     1489176 :                             u               , v                  , q             , dse          ,               &
         141           0 :                             tautmsx         , tautmsy            , dtk           , topflx       , errstring   , &
         142     1489176 :                             tauresx         , tauresy            , itaures       , cpairv       , dse_top, &
         143             :                             do_molec_diff   , use_temperature_molec_diff, vd_lu_qdecomp, &
         144     5956704 :                             ubc_mmr, ubc_flux, kvt, pmid, &
         145     4467528 :                             cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx, nbot_molec, &
         146     1489176 :                             kq_scal, mw_fac)
         147             :
         148             :     !-------------------------------------------------------------------------- !
         149             :     ! Driver routine to compute vertical diffusion of momentum, moisture, trace !
         150             :     ! constituents and dry static energy. The new temperature is computed from  !
         151             :     ! the diffused dry static energy.                                           !
         152             :     ! Turbulent diffusivities and boundary layer nonlocal transport terms are   !
         153             :     ! obtained from the turbulence module.                                      !
         154             :     !-------------------------------------------------------------------------- !
         155             :
         156             : !    Used for CAM debugging.
         157             : !    use phys_debug_util,    only : phys_debug_col
         158             : !    use time_manager,       only : is_first_step, get_nstep
         159             :
         160             :     use coords_1d, only: Coords1D
         161             :     use linear_1d_operators, only : BoundaryType, BoundaryFixedLayer, &
         162             :          BoundaryData, BoundaryFlux, TriDiagDecomp
         163             :     use vdiff_lu_solver,     only : fin_vol_lu_decomp
         164             :     use vertical_diffusion_solver, only : fin_vol_solve
         165             :     use beljaars_drag_cam,   only : do_beljaars
         166             :     ! FIXME: This should not be needed
         167             :     use air_composition,     only: rairv
         168             :
         169             :     use phys_control,        only : phys_getopts
         170             :
         171             :   ! Modification : Ideally, we should diffuse 'liquid-ice static energy' (sl), not the dry static energy.
         172             :   !                Also, vertical diffusion of cloud droplet number concentration and aerosol number
         173             :   !                concentration should be done very carefully in the future version.
         174             :
         175             :     ! --------------- !
         176             :     ! Input Arguments !
         177             :     ! --------------- !
         178             :
         179             :     integer,  intent(in)    :: lchnk
         180             :     integer,  intent(in)    :: pcols
         181             :     integer,  intent(in)    :: pver
         182             :     integer,  intent(in)    :: ncnst
         183             :     integer,  intent(in)    :: ncol                      ! Number of atmospheric columns
         184             :     integer,  intent(in)    :: itaures                   ! Indicator determining whether 'tauresx,tauresy'
         185             :                                                          ! is updated (1) or non-updated (0) in this subroutine.
         186             :
         187             :     type(Coords1D), intent(in) :: p                      ! Pressure coordinates [ Pa ]
         188             :     real(r8), intent(in)    :: tint(pcols,pver+1)        ! Temperature [ K ]
         189             :     real(r8), intent(in)    :: t(pcols,pver)             ! Temperature [ K ]
         190             :     real(r8), intent(in)    :: rhoi(pcols,pver+1)        ! Density of air at interfaces [ kg/m3 ]
         191             :     real(r8), intent(in)    :: ztodt                     ! 2 delta-t [ s ]
         192             :     real(r8), intent(in)    :: taux(pcols)               ! Surface zonal      stress.
         193             :                                                          ! Input u-momentum per unit time per unit area into the atmosphere [ N/m2 ]
         194             :     real(r8), intent(in)    :: tauy(pcols)               ! Surface meridional stress.
         195             :                                                          ! Input v-momentum per unit time per unit area into the atmosphere [ N/m2 ]
         196             :     real(r8), intent(in)    :: shflx(pcols)              ! Surface sensible heat flux [ W/m2 ]
         197             :     real(r8), intent(in)    :: cflx(pcols,ncnst)         ! Surface constituent flux [ kg/m2/s ]
         198             :     real(r8), intent(in)    :: zi(pcols,pver+1)          ! Interface heights [ m ]
         199             :     real(r8), intent(in)    :: ksrftms(pcols)            ! Surface drag coefficient for turbulent mountain stress. > 0. [ kg/s/m2 ]
         200             :     real(r8), intent(in)    :: dragblj(pcols,pver)       ! Drag profile from Beljaars SGO form drag  > 0. [ 1/s ]
         201             :     real(r8), intent(in)    :: qmincg(ncnst)             ! Minimum constituent mixing ratios from cg fluxes
         202             :     real(r8), intent(in)    :: cpairv(pcols,pver)        ! Specific heat at constant pressure
         203             :     real(r8), intent(in)    :: kvh(pcols,pver+1)         ! Eddy diffusivity for heat [ m2/s ]
         204             :
         205             :     logical,  intent(in)    :: do_molec_diff             ! Flag indicating multiple constituent diffusivities
         206             :     logical,  intent(in)    :: use_temperature_molec_diff! Flag indicating that molecular diffusion should apply to temperature, not
         207             :                                                          ! dry static energy.
         208             :
         209             :     type(vdiff_selector), intent(in) :: fieldlist        ! Array of flags selecting which fields to diffuse
         210             :     type(vdiff_selector), intent(in) :: fieldlistm       ! Array of flags selecting which fields for molecular diffusion
         211             :
         212             :     ! Dry static energy top boundary condition.
         213             :     real(r8), intent(in) :: dse_top(pcols)
         214             :
         215             :     real(r8), intent(in) :: kvm(pcols,pver+1)         ! Eddy viscosity ( Eddy diffusivity for momentum ) [ m2/s ]
         216             :     real(r8), intent(in) :: kvq(pcols,pver+1)         ! Eddy diffusivity for constituents
         217             :     real(r8), intent(in) :: cgs(pcols,pver+1)         ! Counter-gradient star [ cg/flux ]
         218             :     real(r8), intent(in) :: cgh(pcols,pver+1)         ! Counter-gradient term for heat
         219             :
         220             :     ! ---------------------- !
         221             :     ! Input-Output Arguments !
         222             :     ! ---------------------- !
         223             :
         224             :     real(r8), intent(inout) :: u(pcols,pver)             ! U wind. This input is the 'raw' input wind to
         225             :                                                          ! PBL scheme without iterative provisional update. [ m/s ]
         226             :     real(r8), intent(inout) :: v(pcols,pver)             ! V wind. This input is the 'raw' input wind to PBL scheme
         227             :                                                          ! without iterative provisional update. [ m/s ]
         228             :     real(r8), intent(inout) :: q(pcols,pver,ncnst)       ! Moisture and trace constituents [ kg/kg, #/kg ? ]
         229             :     real(r8), intent(inout) :: dse(pcols,pver)           ! Dry static energy [ J/kg ]
         230             :
         231             :     real(r8), intent(inout) :: tauresx(pcols)            ! Input  : Reserved surface stress at previous time step
         232             :     real(r8), intent(inout) :: tauresy(pcols)            ! Output : Reserved surface stress at current  time step
         233             :
         234             :     ! ---------------- !
         235             :     ! Output Arguments !
         236             :     ! ---------------- !
         237             :
         238             :     real(r8), intent(out)   :: dtk(pcols,pver)           ! T tendency from KE dissipation
         239             :     real(r8), intent(out)   :: tautmsx(pcols)            ! Implicit zonal      turbulent mountain surface stress
         240             :                                                          ! [ N/m2 = kg m/s /s/m2 ]
         241             :     real(r8), intent(out)   :: tautmsy(pcols)            ! Implicit meridional turbulent mountain surface stress
         242             :                                                          ! [ N/m2 = kg m/s /s/m2 ]
         243             :     real(r8), intent(out)   :: topflx(pcols)             ! Molecular heat flux at the top interface
         244             :     character(128), intent(out) :: errstring             ! Output status
         245             :
         246             :     ! ------------------ !
         247             :     ! Optional Arguments !
         248             :     ! ------------------ !
         249             :
         250             :     ! The molecular diffusion module will likely change significantly in
         251             :     ! the future, and this module may directly depend on it after that.
         252             :     ! Until then, we have these highly specific interfaces hard-coded.
         253             :
         254             :     optional :: vd_lu_qdecomp        ! Constituent-dependent molecular diffusivity routine
         255             :
         256             :     interface
         257             :        function vd_lu_qdecomp( &
         258             :             pcols , pver   , ncol       , fixed_ubc  , mw     , &
         259             :             kv    , kq_scal, mw_facm    , dpidz_sq   , coords , &
         260             :             interface_boundary, molec_boundary, &
         261             :             tint  , ztodt  , nbot_molec , &
         262             :             lchnk , t          , m      , no_molec_decomp) result(decomp)
         263             :          import
         264             :          integer,  intent(in)    :: pcols
         265             :          integer,  intent(in)    :: pver
         266             :          integer,  intent(in)    :: ncol
         267             :          integer,  intent(in)    :: nbot_molec
         268             :          logical,  intent(in)    :: fixed_ubc
         269             :          real(r8), intent(in)    :: kv(pcols,pver+1)
         270             :          real(r8), intent(in)    :: kq_scal(pcols,pver+1)
         271             :          real(r8), intent(in)    :: mw
         272             :          real(r8), intent(in)    :: mw_facm(pcols,pver+1)
         273             :          real(r8), intent(in)    :: dpidz_sq(ncol,pver+1)
         274             :          type(Coords1D), intent(in) :: coords
         275             :          type(BoundaryType), intent(in) :: interface_boundary
         276             :          type(BoundaryType), intent(in) :: molec_boundary
         277             :          real(r8), intent(in)    :: tint(pcols,pver+1)
         278             :          real(r8), intent(in)    :: ztodt
         279             :          integer,  intent(in)    :: lchnk
         280             :          real(r8), intent(in)    :: t(pcols,pver)
         281             :          integer,  intent(in)    :: m
         282             :          type(TriDiagDecomp), intent(in) :: no_molec_decomp
         283             :          type(TriDiagDecomp) :: decomp
         284             :        end function vd_lu_qdecomp
         285             :     end interface
         286             :
         287             :     real(r8), intent(in), optional :: ubc_mmr(pcols,ncnst) ! Upper boundary mixing ratios [ kg/kg ]
         288             :     real(r8), intent(in), optional :: ubc_flux(pcols,ncnst)      ! Upper boundary flux [ kg/s/m^2 ]
         289             :
         290             :     real(r8), intent(in), optional :: kvt(pcols,pver+1) ! Kinematic molecular conductivity
         291             :
         292             :     ! FIXME: This input should not be needed (and should not be passed in in vertical_diffusion).
         293             :     real(r8), intent(in), optional :: pmid(pcols,pver)
         294             :
         295             :     real(r8), intent(in), optional :: cnst_mw(ncnst)          ! Molecular weight [ kg/kmole ]
         296             :     logical, intent(in), optional  :: cnst_fixed_ubc(ncnst)   ! Whether upper boundary condition is fixed
         297             :     logical, intent(in), optional  :: cnst_fixed_ubflx(ncnst) ! Whether upper boundary flux is a fixed non-zero value
         298             :
         299             :     integer, intent(in), optional  :: nbot_molec              ! Bottom level where molecular diffusivity is applied
         300             :
         301             :     ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity
         302             :     real(r8), intent(in), optional :: kq_scal(pcols,pver+1)
         303             :     ! Local sqrt(1/M_q + 1/M_d) for each constituent
         304             :     real(r8), intent(in), optional :: mw_fac(pcols,pver+1,ncnst)
         305             :
         306             :     ! --------------- !
         307             :     ! Local Variables !
         308             :     ! --------------- !
         309             :
         310             :     integer  :: i, k, m                                  ! Longitude, level, constituent indices
         311     2978352 :     logical  :: lqtst(pcols)                             ! Adjust vertical profiles
         312             :
         313             :     ! LU decomposition information.
         314     1489176 :     type(TriDiagDecomp) :: decomp
         315     1489176 :     type(TriDiagDecomp) :: no_molec_decomp
         316             :
         317             :     ! Square of derivative of pressure with height (on interfaces).
         318     2978352 :     real(r8) :: dpidz_sq(ncol,pver+1)
         319             :
         320             :     ! Pressure coordinates over the molecular diffusion range only.
         321     1489176 :     type(Coords1D) :: p_molec
         322             :
         323             :     ! Boundary layer objects
         324     1489176 :     type(BoundaryType) :: interface_boundary
         325     1489176 :     type(BoundaryType) :: molec_boundary
         326             :
         327     2978352 :     real(r8) :: tmp1(pcols)                              ! Temporary storage
         328     2978352 :     real(r8) :: tmpi1(pcols,pver+1)                      ! Interface KE dissipation
         329     2978352 :     real(r8) :: tmpi2(pcols,pver+1)                      ! dt*(g*rho)**2/dp at interfaces
         330     2978352 :     real(r8) :: keg_in(pcols,pver)                       ! KE on entry to subroutine
         331     2978352 :     real(r8) :: keg_out(pcols,pver)                      ! KE after U and V dissipation/diffusion
         332     2978352 :     real(r8) :: rrho(pcols)                              ! 1./bottom level density
         333             :
         334     2978352 :     real(r8) :: tautotx(pcols)                           ! Total surface stress ( zonal )
         335     2978352 :     real(r8) :: tautoty(pcols)                           ! Total surface stress ( meridional )
         336             :
         337     2978352 :     real(r8) :: dinp_u(pcols,pver+1)                     ! Vertical difference at interfaces, input u
         338     2978352 :     real(r8) :: dinp_v(pcols,pver+1)                     ! Vertical difference at interfaces, input v
         339             :     real(r8) :: dout_u                                   ! Vertical difference at interfaces, output u
         340             :     real(r8) :: dout_v                                   ! Vertical difference at interfaces, output v
         341             :
         342     2978352 :     real(r8) :: qtm(pcols,pver)                          ! Temporary copy of q
         343             :
         344     2978352 :     real(r8) :: ws(pcols)                                ! Lowest-level wind speed [ m/s ]
         345     2978352 :     real(r8) :: tau(pcols)                               ! Turbulent surface stress ( not including mountain stress )
         346     2978352 :     real(r8) :: ksrfturb(pcols)                          ! Surface drag coefficient of 'normal' stress. > 0.
         347             :                                                          ! Virtual mass input per unit time per unit area [ kg/s/m2 ]
         348     2978352 :     real(r8) :: ksrf(pcols)                              ! Surface drag coefficient of 'normal' stress +
         349             :                                                          ! Surface drag coefficient of 'tms' stress.  > 0. [ kg/s/m2 ]
         350     2978352 :     real(r8) :: usum_in(pcols)                           ! Vertical integral of input u-momentum. Total zonal
         351             :                                                          ! momentum per unit area in column  [ sum of u*dp/g = kg m/s m-2 ]
         352     2978352 :     real(r8) :: vsum_in(pcols)                           ! Vertical integral of input v-momentum. Total meridional
         353             :                                                          ! momentum per unit area in column [ sum of v*dp/g = kg m/s m-2 ]
         354     2978352 :     real(r8) :: usum_mid(pcols)                          ! Vertical integral of u-momentum after adding explicit residual stress
         355     2978352 :     real(r8) :: vsum_mid(pcols)                          ! Vertical integral of v-momentum after adding explicit residual stress
         356     2978352 :     real(r8) :: usum_out(pcols)                          ! Vertical integral of u-momentum after doing implicit diffusion
         357     2978352 :     real(r8) :: vsum_out(pcols)                          ! Vertical integral of v-momentum after doing implicit diffusion
         358     2978352 :     real(r8) :: tauimpx(pcols)                           ! Actual net stress added at the current step other than mountain stress
         359     2978352 :     real(r8) :: tauimpy(pcols)                           ! Actual net stress added at the current step other than mountain stress
         360             :     real(r8) :: ramda                                    ! dt/timeres [ no unit ]
         361             :
         362     2978352 :     real(r8) :: taubljx(pcols)                           ! recomputed explicit/residual beljaars stress
         363     2978352 :     real(r8) :: taubljy(pcols)                           ! recomputed explicit/residual beljaars stress
         364             :
         365             :     ! Rate at which external (surface) stress damps wind speeds (1/s).
         366     2978352 :     real(r8) :: tau_damp_rate(ncol, pver)
         367             :
         368             :     ! Combined molecular and eddy diffusion.
         369     2978352 :     real(r8) :: kv_total(pcols,pver+1)
         370             :
         371             :     logical  :: use_spcam
         372             :
         373             :     !--------------------------------
         374             :     ! Variables needed for WACCM-X
         375             :     !--------------------------------
         376     2978352 :     real(r8) :: ttemp(ncol,pver)                         ! temporary temperature array
         377     1489176 :     real(r8) :: ttemp0(ncol,pver)                        ! temporary temperature array
         378             :
         379             :     ! ------------------------------------------------ !
         380             :     ! Parameters for implicit surface stress treatment !
         381             :     ! ------------------------------------------------ !
         382             :
         383             :     real(r8), parameter :: wsmin = 1._r8                 ! Minimum sfc wind speed for estimating frictional
         384             :                                                          ! transfer velocity ksrf. [ m/s ]
         385             :     real(r8), parameter :: ksrfmin = 1.e-4_r8            ! Minimum surface drag coefficient [ kg/s/m^2 ]
         386             :     real(r8), parameter :: timeres = 7200._r8            ! Relaxation time scale of residual stress ( >= dt ) [ s ]
         387             :
         388             :     ! ----------------------- !
         389             :     ! Main Computation Begins !
         390             :     ! ----------------------- !
         391             :
         392     1489176 :     call phys_getopts(use_spcam_out = use_spcam)
         393             :
         394     1489176 :     errstring = ''
         395     1489176 :     if( ( diffuse(fieldlist,'u') .or. diffuse(fieldlist,'v') ) .and. .not. diffuse(fieldlist,'s') ) then
         396           0 :           errstring = 'diffusion_solver.compute_vdiff: must diffuse s if diffusing u or v'
         397           0 :           return
         398             :     end if
         399             :
         400             :     !--------------------------------------- !
         401             :     ! Computation of Molecular Diffusivities !
         402             :     !--------------------------------------- !
         403             :
         404             :   ! Modification : Why 'kvq' is not changed by molecular diffusion ?
         405             :
         406     1489176 :     if( do_molec_diff ) then
         407             :
         408             :         if( (.not.present(vd_lu_qdecomp)) .or. (.not.present(kvt)) &
         409           0 :              .or. (.not. present(ubc_mmr)) .or. (.not. present(ubc_flux)) ) then
         410           0 :               errstring = 'compute_vdiff: do_molec_diff true but vd_lu_qdecomp or kvt missing'
         411           0 :               return
         412             :         endif
         413             :
         414           0 :         p_molec = p%section([1, ncol], [1, nbot_molec])
         415           0 :         molec_boundary = BoundaryFixedLayer(p%del(:,nbot_molec+1))
         416             :
         417             :     endif
         418             :
         419             :     ! Boundary condition for a fixed concentration directly on a boundary
         420             :     ! interface (i.e. a boundary layer of size 0).
         421     1489176 :     interface_boundary = BoundaryFixedLayer(spread(0._r8, 1, ncol))
         422             :
         423             :     ! Note that the *derivative* dp/dz is g*rho
         424  2338872120 :     dpidz_sq = gravit*rhoi(:ncol,:)
         425  2338872120 :     dpidz_sq = dpidz_sq * dpidz_sq
         426             :
         427    24865776 :     rrho(:ncol) = rair  * t(:ncol,pver) / p%mid(:,pver)
         428             :
         429    24865776 :     tmpi2(:ncol,1) = ztodt * dpidz_sq(:,1) / ( p%mid(:,1) - p%ifc(:,1) )
         430  2289140568 :     tmpi2(:ncol,2:pver) = ztodt * dpidz_sq(:,2:pver) * p%rdst
         431             :
         432             :     ! FIXME: The following four lines are kept in only to preserve answers;
         433             :     !        they really should be taken out completely.
         434     1489176 :     if (do_molec_diff) &
         435           0 :          tmpi2(:ncol,1) = ztodt * (gravit * rhoi(:ncol,1))**2 / ( pmid(:ncol,1) - p%ifc(:,1) )
         436    24865776 :     dpidz_sq(:,1) = gravit*(p%ifc(:,1) / (rairv(:ncol,1,lchnk)*t(:ncol,1)))
         437    24865776 :     dpidz_sq(:,1) = dpidz_sq(:,1)*dpidz_sq(:,1)
         438             :
         439    24865776 :     tmp1(:ncol) = ztodt * gravit * p%rdel(:,pver)
         440             :
         441             :     !---------------------------- !
         442             :     ! Diffuse Horizontal Momentum !
         443             :     !---------------------------- !
         444             :
         445   139982544 :     do k = 1, pver
         446  2314006344 :        do i = 1, ncol
         447  2312517168 :           keg_in(i,k) = 0.5_r8 * ( u(i,k)*u(i,k) + v(i,k)*v(i,k) )
         448             :        end do
         449             :     end do
         450             :
         451     1489176 :     if( diffuse(fieldlist,'u') .or. diffuse(fieldlist,'v') ) then
         452             :
         453             :         ! Compute the vertical upward differences of the input u,v for KE dissipation
         454             :         ! at each interface.
         455             :         ! Velocity = 0 at surface, so difference at the bottom interface is -u,v(pver)
         456             :         ! These 'dinp_u, dinp_v' are computed using the non-diffused input wind.
         457             :
         458    24865776 :         do i = 1, ncol
         459    23376600 :            dinp_u(i,1) = 0._r8
         460    23376600 :            dinp_v(i,1) = 0._r8
         461    23376600 :            dinp_u(i,pver+1) = -u(i,pver)
         462    24865776 :            dinp_v(i,pver+1) = -v(i,pver)
         463             :         end do
         464   138493368 :         do k = 2, pver
         465  2289140568 :            do i = 1, ncol
         466  2150647200 :               dinp_u(i,k) = u(i,k) - u(i,k-1)
         467  2287651392 :               dinp_v(i,k) = v(i,k) - v(i,k-1)
         468             :            end do
         469             :         end do
         470             :
         471             :        ! -------------------------------------------------------------- !
         472             :        ! Do 'Implicit Surface Stress' treatment for numerical stability !
         473             :        ! in the lowest model layer.                                     !
         474             :        ! -------------------------------------------------------------- !
         475             :
         476     1489176 :        if( do_iss ) then
         477             :
         478             :          ! Compute surface drag coefficient for implicit diffusion
         479             :          ! including turbulent mountain stress.
         480             :
         481    24865776 :            do i = 1, ncol
         482    23376600 :               ws(i)       = max( sqrt( u(i,pver)**2._r8 + v(i,pver)**2._r8 ), wsmin )
         483    23376600 :               tau(i)      = sqrt( taux(i)**2._r8 + tauy(i)**2._r8 )
         484    24865776 :               ksrfturb(i) = max( tau(i) / ws(i), ksrfmin )
         485             :            end do
         486    24865776 :            ksrf(:ncol) = ksrfturb(:ncol) + ksrftms(:ncol)  ! Do all surface stress ( normal + tms ) implicitly
         487             :
         488             :          ! Vertical integration of input momentum.
         489             :          ! This is total horizontal momentum per unit area [ kg*m/s/m2 ] in each column.
         490             :          ! Note (u,v) are the raw input to the PBL scheme, not the
         491             :          ! provisionally-marched ones within the iteration loop of the PBL scheme.
         492             :
         493    24865776 :            do i = 1, ncol
         494    23376600 :               usum_in(i) = 0._r8
         495    23376600 :               vsum_in(i) = 0._r8
         496  2198889576 :               do k = 1, pver
         497  2174023800 :                  usum_in(i) = usum_in(i) + (1._r8/gravit)*u(i,k)*p%del(i,k)
         498  2197400400 :                  vsum_in(i) = vsum_in(i) + (1._r8/gravit)*v(i,k)*p%del(i,k)
         499             :               end do
         500             :            end do
         501             :
         502             :          ! Add residual stress of previous time step explicitly into the lowest
         503             :          ! model layer with a relaxation time scale of 'timeres'.
         504             :
         505     1489176 :            if (am_correction) then
         506             :               ! preserve time-mean torque
         507             :               ramda         = 1._r8
         508             :            else
         509     1489176 :               ramda         = ztodt / timeres
         510             :            endif
         511             :
         512    24865776 :            u(:ncol,pver) = u(:ncol,pver) + tmp1(:ncol)*tauresx(:ncol)*ramda
         513    24865776 :            v(:ncol,pver) = v(:ncol,pver) + tmp1(:ncol)*tauresy(:ncol)*ramda
         514             :
         515             :          ! Vertical integration of momentum after adding explicit residual stress
         516             :          ! into the lowest model layer.
         517             :
         518    24865776 :            do i = 1, ncol
         519    23376600 :               usum_mid(i) = 0._r8
         520    23376600 :               vsum_mid(i) = 0._r8
         521  2198889576 :               do k = 1, pver
         522  2174023800 :                  usum_mid(i) = usum_mid(i) + (1._r8/gravit)*u(i,k)*p%del(i,k)
         523  2197400400 :                  vsum_mid(i) = vsum_mid(i) + (1._r8/gravit)*v(i,k)*p%del(i,k)
         524             :               end do
         525             :            end do
         526             :
         527             :        else
         528             :
         529             :          ! In this case, do 'turbulent mountain stress' implicitly,
         530             :          ! but do 'normal turbulent stress' explicitly.
         531             :          ! In this case, there is no 'residual stress' as long as 'tms' is
         532             :          ! treated in a fully implicit way, which is true.
         533             :
         534             :          ! 1. Do 'tms' implicitly
         535             :
         536           0 :            ksrf(:ncol) = ksrftms(:ncol)
         537             :
         538             :          ! 2. Do 'normal stress' explicitly
         539             :
         540           0 :            u(:ncol,pver) = u(:ncol,pver) + tmp1(:ncol)*taux(:ncol)
         541           0 :            v(:ncol,pver) = v(:ncol,pver) + tmp1(:ncol)*tauy(:ncol)
         542             :
         543             :        end if  ! End of 'do iss' ( implicit surface stress )
         544             :
         545             :        ! --------------------------------------------------------------------------------------- !
         546             :        ! Diffuse horizontal momentum implicitly using tri-diagnonal matrix.                      !
         547             :        ! The 'u,v' are input-output: the output 'u,v' are implicitly diffused winds.             !
         548             :        !    For implicit 'normal' stress : ksrf = ksrftms + ksrfturb,                            !
         549             :        !                                   u(pver) : explicitly include 'residual normal' stress !
         550             :        !    For explicit 'normal' stress : ksrf = ksrftms                                        !
         551             :        !                                   u(pver) : explicitly include 'normal' stress          !
         552             :        ! Note that in all the two cases above, 'tms' is fully implicitly treated.                !
         553             :        ! --------------------------------------------------------------------------------------- !
         554             :
         555             :        ! In most layers, no damping at all.
         556  2314006344 :        tau_damp_rate = 0._r8
         557             :
         558             :        ! Physical interpretation:
         559             :        ! ksrf is stress per unit wind speed.
         560             :        ! p%del / gravit is approximately the mass in the layer per unit of
         561             :        ! surface area.
         562             :        ! Therefore, gravit*ksrf/p%del is the acceleration of wind per unit
         563             :        ! wind speed, i.e. the rate at which wind is exponentially damped by
         564             :        ! surface stress.
         565             :
         566             :        ! Beljaars et al SGO scheme incorporated here. It
         567             :        ! appears as a "3D" tau_damp_rate specification.
         568             :
         569    24865776 :        tau_damp_rate(:,pver) = -gravit*ksrf(:ncol)*p%rdel(:,pver)
         570   139982544 :        do k=1,pver
         571  2314006344 :           tau_damp_rate(:,k) = tau_damp_rate(:,k) + dragblj(:ncol,k)
         572             :        end do
         573             :
         574             :        v(:ncol,:) = fin_vol_solve(ztodt, p, v(:ncol,:), ncol, pver, &
         575             :                          coef_q=tau_damp_rate,                      &
         576  5193553248 :                          coef_q_diff=kvm(:ncol,:)*dpidz_sq)
         577             :
         578             :        u(:ncol,:) = fin_vol_solve(ztodt, p, u(:ncol,:), ncol, pver, &
         579             :                          coef_q=tau_damp_rate,                      &
         580  5193553248 :                          coef_q_diff=kvm(:ncol,:)*dpidz_sq)
         581             :
         582             :
         583             :
         584             :        ! ---------------------------------------------------------------------- !
         585             :        ! Calculate 'total' ( tautotx ) and 'tms' ( tautmsx ) stresses that      !
         586             :        ! have been actually added into the atmosphere at the current time step. !
         587             :        ! Also, update residual stress, if required.                             !
         588             :        ! ---------------------------------------------------------------------- !
         589             :
         590    24865776 :        do i = 1, ncol
         591             :
         592             :           ! Compute the implicit 'tms' using the updated winds.
         593             :           ! Below 'tautmsx(i),tautmsy(i)' are pure implicit mountain stresses
         594             :           ! that has been actually added into the atmosphere both for explicit
         595             :           ! and implicit approach.
         596             :
         597    23376600 :           tautmsx(i) = -ksrftms(i)*u(i,pver)
         598    23376600 :           tautmsy(i) = -ksrftms(i)*v(i,pver)
         599             :
         600             :           ! We want to add vertically-integrated Beljaars drag to residual stress.
         601             :           ! So this has to be calculated locally.
         602             :           ! We may want to rethink the residual drag calculation performed here on. (jtb)
         603    23376600 :           taubljx(i) = 0._r8
         604    23376600 :           taubljy(i) = 0._r8
         605  2197400400 :           do k = 1, pver
         606  2174023800 :              taubljx(i) = taubljx(i) + (1._r8/gravit)*dragblj(i,k)*u(i,k)*p%del(i,k)
         607  2197400400 :              taubljy(i) = taubljy(i) + (1._r8/gravit)*dragblj(i,k)*v(i,k)*p%del(i,k)
         608             :           end do
         609             :
         610    24865776 :           if( do_iss ) then
         611             :
         612             :             ! Compute vertical integration of final horizontal momentum
         613             :
         614    23376600 :               usum_out(i) = 0._r8
         615    23376600 :               vsum_out(i) = 0._r8
         616  2197400400 :               do k = 1, pver
         617  2174023800 :                  usum_out(i) = usum_out(i) + (1._r8/gravit)*u(i,k)*p%del(i,k)
         618  2197400400 :                  vsum_out(i) = vsum_out(i) + (1._r8/gravit)*v(i,k)*p%del(i,k)
         619             :               end do
         620             :
         621             :             ! Compute net stress added into the atmosphere at the current time step.
         622             :             ! Note that the difference between 'usum_in' and 'usum_out' are induced
         623             :             ! by 'explicit residual stress + implicit total stress' for implicit case, while
         624             :             ! by 'explicit normal   stress + implicit tms   stress' for explicit case.
         625             :             ! Here, 'tautotx(i)' is net stress added into the air at the current time step.
         626             :
         627    23376600 :               tauimpx(i) = ( usum_out(i) - usum_in(i) ) / ztodt
         628    23376600 :               tauimpy(i) = ( vsum_out(i) - vsum_in(i) ) / ztodt
         629             :
         630    23376600 :               tautotx(i) = tauimpx(i)
         631    23376600 :               tautoty(i) = tauimpy(i)
         632             :
         633             :             ! Compute residual stress and update if required.
         634             :             ! Note that the total stress we should have added at the current step is
         635             :             ! the sum of 'taux(i) - ksrftms(i)*u(i,pver) + tauresx(i)'.
         636             :
         637    23376600 :               if( itaures .eq. 1 ) then
         638    23376600 :                  tauresx(i) = taux(i) + tautmsx(i) + taubljx(i) + tauresx(i)- tauimpx(i)
         639    23376600 :                  tauresy(i) = tauy(i) + tautmsy(i) + taubljy(i) + tauresy(i)- tauimpy(i)
         640             :               endif
         641             :
         642             :           else
         643             :
         644           0 :              tautotx(i) = tautmsx(i) + taux(i)
         645           0 :              tautoty(i) = tautmsy(i) + tauy(i)
         646           0 :              tauresx(i) = 0._r8
         647           0 :              tauresy(i) = 0._r8
         648             :
         649             :           end if  ! End of 'do_iss' if
         650             :
         651             :        end do ! End of 'do i = 1, ncol' loop
         652             :
         653             :        ! ------------------------------------ !
         654             :        ! Calculate kinetic energy dissipation !
         655             :        ! ------------------------------------ !
         656             :
         657             :      ! Modification : In future, this should be set exactly same as
         658             :      !                the ones in the convection schemes
         659             :
         660             :        ! 1. Compute dissipation term at interfaces
         661             :        !    Note that 'u,v' are already diffused wind, and 'tautotx,tautoty' are
         662             :        !    implicit stress that has been actually added. On the other hand,
         663             :        !    'dinp_u, dinp_v' were computed using non-diffused input wind.
         664             :
         665             :      ! Modification : I should check whether non-consistency between 'u' and 'dinp_u'
         666             :      !                is correctly intended approach. I think so.
         667             :
         668    24865776 :        k = pver + 1
         669    24865776 :        do i = 1, ncol
         670    23376600 :           tmpi1(i,1) = 0._r8
         671             :           tmpi1(i,k) = 0.5_r8 * ztodt * gravit * &
         672    24865776 :                        ( (-u(i,k-1) + dinp_u(i,k))*tautotx(i) + (-v(i,k-1) + dinp_v(i,k))*tautoty(i) )
         673             :        end do
         674             :
         675   138493368 :        do k = 2, pver
         676  2289140568 :           do i = 1, ncol
         677  2150647200 :              dout_u = u(i,k) - u(i,k-1)
         678  2150647200 :              dout_v = v(i,k) - v(i,k-1)
         679  2150647200 :              tmpi1(i,k) = 0.25_r8 * tmpi2(i,k) * kvm(i,k) * &
         680  2287651392 :                           ( dout_u**2 + dout_v**2 + dout_u*dinp_u(i,k) + dout_v*dinp_v(i,k) )
         681             :           end do
         682             :        end do
         683             :
         684     1489176 :        if (do_beljaars) then
         685             :
         686             :           ! 2. Add Kinetic Energy change across dissipation to Static Energy
         687   139982544 :           do k = 1, pver
         688  2314006344 :              do i = 1, ncol
         689  2312517168 :                 keg_out(i,k) = 0.5_r8 * ( u(i,k)*u(i,k) + v(i,k)*v(i,k) )
         690             :              end do
         691             :           end do
         692             :
         693   139982544 :           do k = 1, pver
         694  2314006344 :              do i = 1, ncol
         695  2174023800 :                 dtk(i,k) = keg_in(i,k) - keg_out(i,k)
         696  2312517168 :                 dse(i,k) = dse(i,k) + dtk(i,k) ! + dkeblj(i,k)
         697             :              end do
         698             :           end do
         699             :
         700             :        else
         701             :
         702             :           ! 2. Compute dissipation term at midpoints, add to dry static energy
         703           0 :           do k = 1, pver
         704           0 :              do i = 1, ncol
         705           0 :                 dtk(i,k) = ( tmpi1(i,k+1) + tmpi1(i,k) ) * p%rdel(i,k)
         706           0 :                 dse(i,k) = dse(i,k) + dtk(i,k)
         707             :              end do
         708             :           end do
         709             :
         710             :        end if
         711             :
         712             :     end if ! End of diffuse horizontal momentum, diffuse(fieldlist,'u') routine
         713             :
         714             :     !-------------------------- !
         715             :     ! Diffuse Dry Static Energy !
         716             :     !-------------------------- !
         717             :
         718             :   ! Modification : In future, we should diffuse the fully conservative
         719             :   !                moist static energy,not the dry static energy.
         720             :
         721     1489176 :     if( diffuse(fieldlist,'s') ) then
         722     1489176 :       if (.not. use_spcam) then
         723             :
         724             :        ! Add counter-gradient to input static energy profiles
         725             :
         726   139982544 :          do k = 1, pver
         727   276986736 :             dse(:ncol,k) = dse(:ncol,k) + ztodt * p%rdel(:,k) * gravit  *                &
         728   138493368 :                                         ( rhoi(:ncol,k+1) * kvh(:ncol,k+1) * cgh(:ncol,k+1) &
         729  2729486448 :                                         - rhoi(:ncol,k  ) * kvh(:ncol,k  ) * cgh(:ncol,k  ) )
         730             :          end do
         731             :        endif
         732             :        ! Add the explicit surface fluxes to the lowest layer
         733    24865776 :        dse(:ncol,pver) = dse(:ncol,pver) + tmp1(:ncol) * shflx(:ncol)
         734             :
         735             :      ! Diffuse dry static energy
         736             :
         737             :        !---------------------------------------------------
         738             :        ! Solve for temperature using thermal conductivity
         739             :        !---------------------------------------------------
         740     1489176 :        if ( use_temperature_molec_diff ) then
         741             :           !----------------------------------------------------------------------------------------------------
         742             :           ! In Extended WACCM, kvt is calculated rather kvh. This is because molecular diffusion operates on
         743             :           ! temperature, while eddy diffusion operates on dse.  Also, pass in constituent dependent "constants"
         744             :           !----------------------------------------------------------------------------------------------------
         745             :
         746             :           ! Boundary layer thickness of "0._r8" signifies that the boundary
         747             :           ! condition is defined directly on the top interface.
         748             :
         749           0 :           if (.not. use_spcam) then
         750             :              dse(:ncol,:) = fin_vol_solve(ztodt, p, dse(:ncol,:), ncol, pver, &
         751             :                                   coef_q_diff=kvh(:ncol,:)*dpidz_sq,          &
         752             :                                   upper_bndry=interface_boundary,             &
         753           0 :                                   l_cond=BoundaryData(dse_top(:ncol)))
         754             :           endif
         755             :
         756             :           ! Calculate flux at top interface
         757             :
         758             :           ! Modification : Why molecular diffusion does not work for dry static energy in all layers ?
         759             :
         760             :           topflx(:ncol) =  - kvh(:ncol,1) * tmpi2(:ncol,1) / (ztodt*gravit) * &
         761           0 :                ( dse(:ncol,1) - dse_top(:ncol) )
         762             :
         763           0 :           ttemp0 = t(:ncol,:)
         764           0 :           ttemp = ttemp0
         765             :
         766             :           ! upper boundary is zero flux for extended model
         767           0 :           if (.not. use_spcam) then
         768             :              ttemp = fin_vol_solve(ztodt, p, ttemp, ncol, pver, &
         769           0 :                           coef_q_diff=kvt(:ncol,:)*dpidz_sq,    &
         770           0 :                           coef_q_weight=cpairv(:ncol,:))
         771             :           end if
         772             :
         773             :
         774             :           !-------------------------------------
         775             :           !  Update dry static energy
         776             :           !-------------------------------------
         777           0 :           do k = 1,pver
         778           0 :              dse(:ncol,k) = dse(:ncol,k) + &
         779           0 :                   cpairv(:ncol,k)*(ttemp(:,k) - ttemp0(:,k))
         780             :           enddo
         781             :
         782             :        else
         783             :
         784     1489176 :           if (do_molec_diff) then
         785           0 :              kv_total(:ncol,:) = kvh(:ncol,:) + kvt(:ncol,:)/cpair
         786             :           else
         787  2338872120 :              kv_total(:ncol,:) = kvh(:ncol,:)
         788             :           end if
         789             :
         790             :           ! Boundary layer thickness of "0._r8" signifies that the boundary
         791             :           ! condition is defined directly on the top interface.
         792     1489176 :           if (.not. use_spcam) then
         793             :              dse(:ncol,:) = fin_vol_solve(ztodt, p, dse(:ncol,:), ncol, pver, &
         794             :                                  coef_q_diff=kv_total(:ncol,:)*dpidz_sq,      &
         795             :                                  upper_bndry=interface_boundary,              &
         796  5193553248 :                                  l_cond=BoundaryData(dse_top(:ncol)))
         797             :           end if
         798             :
         799             :           ! Calculate flux at top interface
         800             :
         801             :           ! Modification : Why molecular diffusion does not work for dry static energy in all layers ?
         802             :
         803     1489176 :           if( do_molec_diff ) then
         804             :              topflx(:ncol) =  - kv_total(:ncol,1) * tmpi2(:ncol,1) / (ztodt*gravit) * &
         805           0 :                   ( dse(:ncol,1) - dse_top(:ncol) )
         806             :           else
         807    24865776 :              topflx(:ncol) = 0._r8
         808             :           end if
         809             :
         810             :        endif
         811             :
         812             :     endif
         813             :
         814             :     !---------------------------- !
         815             :     ! Diffuse Water Vapor Tracers !
         816             :     !---------------------------- !
         817             :
         818             :   ! Modification : For aerosols, I need to use separate treatment
         819             :   !                for aerosol mass and aerosol number.
         820             :
         821             :     ! Loop through constituents
         822             :
         823             :     no_molec_decomp = fin_vol_lu_decomp(ztodt, p, &
         824  2338872120 :          coef_q_diff=kvq(:ncol,:)*dpidz_sq)
         825             :
         826    17870112 :     do m = 1, ncnst
         827             :
         828    17870112 :        if( diffuse(fieldlist,'q',m) ) then
         829    16380936 :            if (.not. use_spcam) then
         830             :
         831             :               ! Add the nonlocal transport terms to constituents in the PBL.
         832             :               ! Check for neg q's in each constituent and put the original vertical
         833             :               ! profile back if a neg value is found. A neg value implies that the
         834             :               ! quasi-equilibrium conditions assumed for the countergradient term are
         835             :               ! strongly violated.
         836             :
         837 25454069784 :               qtm(:ncol,:pver) = q(:ncol,:pver,m)
         838             :
         839  1539807984 :               do k = 1, pver
         840  1523427048 :                  q(:ncol,k,m) = q(:ncol,k,m) + &
         841  1523427048 :                                 ztodt * p%rdel(:,k) * gravit  * ( cflx(:ncol,m) * rrho(:ncol) ) * &
         842  1523427048 :                               ( rhoi(:ncol,k+1) * kvh(:ncol,k+1) * cgs(:ncol,k+1)                 &
         843 30024350928 :                               - rhoi(:ncol,k  ) * kvh(:ncol,k  ) * cgs(:ncol,k  ) )
         844             :               end do
         845 25454069784 :               lqtst(:ncol) = all(q(:ncol,1:pver,m) >= qmincg(m), 2)
         846  1539807984 :               do k = 1, pver
         847 25454069784 :                  q(:ncol,k,m) = merge( q(:ncol,k,m), qtm(:ncol,k), lqtst(:ncol) )
         848             :               end do
         849             :            endif
         850             :
         851             :            ! Add the explicit surface fluxes to the lowest layer
         852             :
         853   273523536 :            q(:ncol,pver,m) = q(:ncol,pver,m) + tmp1(:ncol) * cflx(:ncol,m)
         854             :
         855             :            ! Diffuse constituents.
         856             :
         857             :            ! This is for solving molecular diffusion of minor species, thus, for WACCM-X, bypass O and O2 (major species)
         858             :            ! Major species diffusion is calculated separately.  -Hanli Liu
         859             :
         860    16380936 :            if( do_molec_diff .and. diffuse(fieldlistm,'q',m)) then
         861             :
         862           0 :               decomp = vd_lu_qdecomp( pcols , pver   , ncol              , cnst_fixed_ubc(m), cnst_mw(m), &
         863             :                    kvq   , kq_scal, mw_fac(:,:,m) ,dpidz_sq          , p_molec, &
         864             :                    interface_boundary, molec_boundary, &
         865             :                    tint  , ztodt  , nbot_molec       , &
         866           0 :                    lchnk , t                , m         , no_molec_decomp)
         867             :
         868             :               ! This to calculate the upper boundary flux of H.    -Hanli Liu
         869           0 :               if ((cnst_fixed_ubflx(m))) then
         870             :
         871             :                  ! ubc_flux is a flux of mass density through space, i.e.:
         872             :                  ! ubc_flux = rho_i * dz/dt = q_i * rho * dz/dt
         873             :                  ! For flux of mmr through pressure level, multiply by g:
         874             :                  ! q_i * rho * gravit * dz/dt = q_i * dp/dt
         875             :
         876             :                  call decomp%left_div(q(:ncol,:,m), &
         877             :                       l_cond=BoundaryFlux( &
         878           0 :                       -gravit*ubc_flux(:ncol,m), ztodt, &
         879           0 :                       p%del(:,1)))
         880             :
         881             :               else
         882             :                  call decomp%left_div(q(:ncol,:,m), &
         883           0 :                       l_cond=BoundaryData(ubc_mmr(:ncol,m)))
         884             :               end if
         885             :
         886           0 :               call decomp%finalize()
         887             :
         888             :            else
         889             :
         890    16380936 :               if (present(cnst_fixed_ubc)) then
         891             :                  ! explicitly set mmr in top layer for cases where molecular diffusion is not active
         892    16380936 :                  if (cnst_fixed_ubc(m)) then
         893           0 :                     q(:ncol,1,m) = ubc_mmr(:ncol,m)
         894             :                  endif
         895             :               end if
         896             :
         897    16380936 :               if (.not. use_spcam) then
         898    16380936 :                  call no_molec_decomp%left_div(q(:ncol,:,m))
         899             :               end if
         900             :
         901             :            end if
         902             :
         903             :        end if
         904             :     end do
         905             :
         906     1489176 :     call no_molec_decomp%finalize()
         907             :
         908    13402584 :   end subroutine compute_vdiff
         909             :
         910             :   ! =============================================================================== !
         911             :   !                                                                                 !
         912             :   ! =============================================================================== !
         913             :
         914       38400 :   character(128) function vdiff_select( fieldlist, name, qindex )
         915             :     ! --------------------------------------------------------------------- !
         916             :     ! This function sets the field with incoming name as one to be diffused !
         917             :     ! --------------------------------------------------------------------- !
         918             :     type(vdiff_selector), intent(inout)        :: fieldlist
         919             :     character(*),         intent(in)           :: name
         920             :     integer,              intent(in), optional :: qindex
         921             :
         922       38400 :     vdiff_select = ''
         923        1536 :     select case (name)
         924             :     case ('u','U')
         925        1536 :        fieldlist%fields(1) = .true.
         926             :     case ('v','V')
         927        1536 :        fieldlist%fields(2) = .true.
         928             :     case ('s','S')
         929        1536 :        fieldlist%fields(3) = .true.
         930             :     case ('q','Q')
         931       33792 :        if( present(qindex) ) then
         932       33792 :            fieldlist%fields(3 + qindex) = .true.
         933             :        else
         934           0 :            fieldlist%fields(4) = .true.
         935             :        endif
         936             :     case default
         937       38400 :        write(vdiff_select,*) 'Bad argument to vdiff_index: ', name
         938             :     end select
         939       38400 :     return
         940             :
         941             :   end function vdiff_select
         942             :
         943           0 :   type(vdiff_selector) function not(a)
         944             :     ! ------------------------------------------------------------- !
         945             :     ! This function extends .not. to operate on type vdiff_selector !
         946             :     ! ------------------------------------------------------------- !
         947             :     type(vdiff_selector), intent(in)  :: a
         948           0 :     allocate(not%fields(size(a%fields)))
         949           0 :     not%fields = .not. a%fields
         950           0 :   end function not
         951             :
         952     2978352 :   logical function my_any(a)
         953             :     ! -------------------------------------------------- !
         954             :     ! This function extends the intrinsic function 'any' !
         955             :     ! to operate on type vdiff_selector                  !
         956             :     ! -------------------------------------------------- !
         957             :     type(vdiff_selector), intent(in) :: a
         958    23826816 :     my_any = any(a%fields)
         959     2978352 :   end function my_any
         960             :
         961    41696928 :   logical function diffuse(fieldlist,name,qindex)
         962             :     ! ---------------------------------------------------------------------------- !
         963             :     ! This function reports whether the field with incoming name is to be diffused !
         964             :     ! ---------------------------------------------------------------------------- !
         965             :     type(vdiff_selector), intent(in)           :: fieldlist
         966             :     character(*),         intent(in)           :: name
         967             :     integer,              intent(in), optional :: qindex
         968             :
         969     2978352 :     select case (name)
         970             :     case ('u','U')
         971     2978352 :        diffuse = fieldlist%fields(1)
         972             :     case ('v','V')
         973     2978352 :        diffuse = fieldlist%fields(2)
         974             :     case ('s','S')
         975     2978352 :        diffuse = fieldlist%fields(3)
         976             :     case ('q','Q')
         977    32761872 :        if( present(qindex) ) then
         978    32761872 :            diffuse = fieldlist%fields(3 + qindex)
         979             :        else
         980           0 :            diffuse = fieldlist%fields(4)
         981             :        endif
         982             :     case default
         983    41696928 :        diffuse = .false.
         984             :     end select
         985             :     return
         986             :   end function diffuse
         987             :
         988           0 : end module diffusion_solver

Generated by: LCOV version 1.14