phys/module_cam_bl_diffusion_solver.F

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