






  module diffusion_solver

  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  



  use module_cam_support,   only : iulog


  implicit none
  private       
  save

  integer, parameter :: r8 = selected_real_kind(12)      

  
  
  

  public init_vdiff                                      
  public compute_vdiff                                   
  public vd_lu_solve                                     
  public vd_lu_decomp                                    
  public vdiff_selector                                  
  public vdiff_select                                    
  public operator(.not.)                                 
  public any                                             

  integer, public   :: nbot_molec                        
 
  

  type vdiff_selector 
       private
       logical, pointer, dimension(:) :: fields
  end type vdiff_selector

  

  interface operator(.not.)
       module procedure not
  end interface

  

  interface any                           
       module procedure my_any
  end interface

  
  
  

  real(r8), private   :: cpair                           
  real(r8), private   :: gravit                          
  real(r8), private   :: rair                            
  real(r8), private   :: zvir                            
  real(r8), private   :: latvap                          
  real(r8), private   :: karman                          

  

  real(r8), parameter :: z0fac   = 0.025_r8              
  real(r8), parameter :: z0max   = 100._r8               
  real(r8), parameter :: horomin = 10._r8                
  real(r8), parameter :: dv2min  = 0.01_r8               
  real(r8), private   :: oroconst                        

  contains

  
  
  

  subroutine init_vdiff( kind, ncnst, rair_in, gravit_in, fieldlist_wet, fieldlist_dry, errstring )

    integer,              intent(in)  :: kind            
    integer,              intent(in)  :: ncnst           
    real(r8),             intent(in)  :: rair_in         
    real(r8),             intent(in)  :: gravit_in       
    type(vdiff_selector), intent(out) :: fieldlist_wet   
    type(vdiff_selector), intent(out) :: fieldlist_dry   
    character(128),       intent(out) :: errstring       
    
    errstring = ''
    if( kind .ne. r8 ) then
        write(iulog,*) 'KIND of reals passed to init_vdiff -- exiting.'

        call wrf_message(iulog)

        errstring = 'init_vdiff'
        return
    endif

    rair   = rair_in     
    gravit = gravit_in 

    allocate( fieldlist_wet%fields( 3 + ncnst ) )
    fieldlist_wet%fields(:) = .false.

    allocate( fieldlist_dry%fields( 3 + ncnst ) )
    fieldlist_dry%fields(:) = .false.

  end subroutine init_vdiff

  
  
  

  subroutine compute_vdiff( lchnk           ,                                                                   &
                            pcols           , pver               , ncnst         , ncol         , pmid        , &
                            pint            , rpdel              , t             , ztodt        , taux        , &
                            tauy            , shflx              , cflx          , ntop         , nbot        , &
                            kvh             , kvm                , kvq           , cgs          , cgh         , &
                            zi              , ksrftms            , qmincg        , fieldlist    ,               &
                            u               , v                  , q             , dse          ,               &
                            tautmsx         , tautmsy            , dtk           , topflx       , errstring   , &
                            tauresx         , tauresy            , itaures       ,                              &
                            do_molec_diff   , compute_molec_diff , vd_lu_qdecomp )

    
    
    
    
    
    
    





  
  
  
  

    
    
    

    integer,  intent(in)    :: lchnk                   
    integer,  intent(in)    :: pcols
    integer,  intent(in)    :: pver
    integer,  intent(in)    :: ncnst
    integer,  intent(in)    :: ncol                      
    integer,  intent(in)    :: ntop                      
    integer,  intent(in)    :: nbot                      
    integer,  intent(in)    :: itaures                   

    real(r8), intent(in)    :: pmid(pcols,pver)          
    real(r8), intent(in)    :: pint(pcols,pver+1)        
    real(r8), intent(in)    :: rpdel(pcols,pver)         
    real(r8), intent(in)    :: t(pcols,pver)             
    real(r8), intent(in)    :: ztodt                     
    real(r8), intent(in)    :: taux(pcols)               
    real(r8), intent(in)    :: tauy(pcols)               
    real(r8), intent(in)    :: shflx(pcols)              
    real(r8), intent(in)    :: cflx(pcols,ncnst)         
    real(r8), intent(in)    :: zi(pcols,pver+1)          
    real(r8), intent(in)    :: ksrftms(pcols)            
    real(r8), intent(in)    :: qmincg(ncnst)             

    logical,  intent(in)         :: do_molec_diff        
    integer,  external, optional :: compute_molec_diff   
    integer,  external, optional :: vd_lu_qdecomp        
    type(vdiff_selector), intent(in) :: fieldlist        

    
    
    

    real(r8), intent(inout) :: kvh(pcols,pver+1)         
    real(r8), intent(inout) :: kvm(pcols,pver+1)         
    real(r8), intent(inout) :: kvq(pcols,pver+1)         
    real(r8), intent(inout) :: cgs(pcols,pver+1)         
    real(r8), intent(inout) :: cgh(pcols,pver+1)         

    real(r8), intent(inout) :: u(pcols,pver)             
    real(r8), intent(inout) :: v(pcols,pver)             
    real(r8), intent(inout) :: q(pcols,pver,ncnst)       
    real(r8), intent(inout) :: dse(pcols,pver)           

    real(r8), intent(inout) :: tauresx(pcols)            
    real(r8), intent(inout) :: tauresy(pcols)            

    
    
    

    real(r8), intent(out)   :: dtk(pcols,pver)           
    real(r8), intent(out)   :: tautmsx(pcols)            
    real(r8), intent(out)   :: tautmsy(pcols)            
    real(r8), intent(out)   :: topflx(pcols)             
    character(128), intent(out) :: errstring             

    
    
    

    integer  :: i, k, m, icol                            
    integer  :: status                                   
    integer  :: ntop_molec                               
    logical  :: lqtst(pcols)                             
    logical  :: need_decomp                              
    logical  :: cnst_fixed_ubc(ncnst)                    
    logical  :: do_iss                                   

    real(r8) :: tmpm(pcols,pver)                         
    real(r8) :: ca(pcols,pver)                           
    real(r8) :: cc(pcols,pver)                           
    real(r8) :: dnom(pcols,pver)                         

    real(r8) :: tmp1(pcols)                              
    real(r8) :: tmpi1(pcols,pver+1)                      
    real(r8) :: tint(pcols,pver+1)                       
    real(r8) :: rhoi(pcols,pver+1)                       
    real(r8) :: tmpi2(pcols,pver+1)                      
    real(r8) :: rrho(pcols)                              

    real(r8) :: zero(pcols)                              
    real(r8) :: tautotx(pcols)                           
    real(r8) :: tautoty(pcols)                           

    real(r8) :: dinp_u(pcols,pver+1)                     
    real(r8) :: dinp_v(pcols,pver+1)                     
    real(r8) :: dout_u                                   
    real(r8) :: dout_v                                   
    real(r8) :: dse_top(pcols)                           
    real(r8) :: cc_top(pcols)                            
    real(r8) :: cd_top(pcols)                            
    real(r8) :: rghd(pcols,pver+1)                       

    real(r8) :: qtm(pcols,pver)                          
    real(r8) :: kq_scal(pcols,pver+1)                    
    real(r8) :: mw_fac(ncnst)                            
    real(r8) :: cnst_mw(ncnst)                           
    real(r8) :: ubc_mmr(pcols,ncnst)                     
    real(r8) :: ubc_t(pcols)                             

    real(r8) :: ws(pcols)                                
    real(r8) :: tau(pcols)                               
    real(r8) :: ksrfturb(pcols)                          
    real(r8) :: ksrf(pcols)                              
    real(r8) :: usum_in(pcols)                           
    real(r8) :: vsum_in(pcols)                           
    real(r8) :: usum_mid(pcols)                          
    real(r8) :: vsum_mid(pcols)                          
    real(r8) :: usum_out(pcols)                          
    real(r8) :: vsum_out(pcols)                          
    real(r8) :: tauimpx(pcols)                           
    real(r8) :: tauimpy(pcols)                           
    real(r8) :: wsmin                                    
    real(r8) :: ksrfmin                                  
    real(r8) :: timeres                                  
    real(r8) :: ramda                                    
    real(r8) :: psum
    real(r8) :: u_in, u_res, tauresx_in
    real(r8) :: v_in, v_res, tauresy_in  

    
    
    

    wsmin    = 1._r8                                     
    ksrfmin  = 1.e-4_r8                                  
    timeres  = 7200._r8                                  



    do_iss = .true. 


    
    
    

    errstring = ''
    if( ( diffuse(fieldlist,'u') .or. diffuse(fieldlist,'v') ) .and. .not. diffuse(fieldlist,'s') ) then
          errstring = 'diffusion_solver.compute_vdiff: must diffuse s if diffusing u or v'
          return
    end if
    zero(:) = 0._r8

    

    tint(:ncol,1) = t(:ncol,1)
    rhoi(:ncol,1) = pint(:ncol,1) / (rair*tint(:ncol,1))
    do k = 2, pver
       do i = 1, ncol
          tint(i,k)  = 0.5_r8 * ( t(i,k) + t(i,k-1) )
          rhoi(i,k)  = pint(i,k) / (rair*tint(i,k))
          tmpi2(i,k) = ztodt * ( gravit*rhoi(i,k) )**2 / ( pmid(i,k) - pmid(i,k-1) )
       end do
    end do
    tint(:ncol,pver+1) = t(:ncol,pver)
    rhoi(:ncol,pver+1) = pint(:ncol,pver+1) / ( rair*tint(:ncol,pver+1) )

    rrho(:ncol) = rair  * t(:ncol,pver) / pmid(:ncol,pver)
    tmp1(:ncol) = ztodt * gravit * rpdel(:ncol,pver)

    
    
    

  

    if( do_molec_diff ) then

        if( (.not.present(compute_molec_diff)) .or. (.not.present(vd_lu_qdecomp)) ) then
              errstring = 'compute_vdiff: do_molec_diff true but compute_molec_diff or vd_lu_qdecomp missing'
              return
        endif

      
      
      
      

        status = compute_molec_diff( lchnk          ,                                                                &
                                     pcols          , pver    , ncnst      , ncol      , t      , pmid   , pint    , &
                                     zi             , ztodt   , kvh        , kvm       , tint   , rhoi   , tmpi2   , &
                                     kq_scal        , ubc_t   , ubc_mmr    , dse_top   , cc_top , cd_top , cnst_mw , &
                                     cnst_fixed_ubc , mw_fac  , ntop_molec , nbot_molec )

    else

        kq_scal(:,:) = 0._r8
        cd_top(:)    = 0._r8
        cc_top(:)    = 0._r8

    endif

    
    
    

    if( diffuse(fieldlist,'u') .or. diffuse(fieldlist,'v') ) then

        
        
        
        

        do i = 1, ncol
           dinp_u(i,1) = 0._r8
           dinp_v(i,1) = 0._r8
           dinp_u(i,pver+1) = -u(i,pver)
           dinp_v(i,pver+1) = -v(i,pver)
        end do
        do k = 2, pver
           do i = 1, ncol
              dinp_u(i,k) = u(i,k) - u(i,k-1)
              dinp_v(i,k) = v(i,k) - v(i,k-1)
           end do
        end do

       
       
       
       

       if( do_iss ) then

         
         

           do i = 1, ncol
              ws(i)       = max( sqrt( u(i,pver)**2._r8 + v(i,pver)**2._r8 ), wsmin )
              tau(i)      = sqrt( taux(i)**2._r8 + tauy(i)**2._r8 )
              ksrfturb(i) = max( tau(i) / ws(i), ksrfmin )
           end do              
           ksrf(:ncol) = ksrfturb(:ncol) + ksrftms(:ncol)  

         
         
         
         

           do i = 1, ncol
              usum_in(i) = 0._r8
              vsum_in(i) = 0._r8
              do k = 1, pver
                 usum_in(i) = usum_in(i) + (1._r8/gravit)*u(i,k)/rpdel(i,k)
                 vsum_in(i) = vsum_in(i) + (1._r8/gravit)*v(i,k)/rpdel(i,k)
              end do
           end do              

         
         

           ramda         = ztodt / timeres
           u(:ncol,pver) = u(:ncol,pver) + tmp1(:ncol)*tauresx(:ncol)*ramda
           v(:ncol,pver) = v(:ncol,pver) + tmp1(:ncol)*tauresy(:ncol)*ramda

         
         

           do i = 1, ncol
              usum_mid(i) = 0._r8
              vsum_mid(i) = 0._r8
              do k = 1, pver
                 usum_mid(i) = usum_mid(i) + (1._r8/gravit)*u(i,k)/rpdel(i,k)
                 vsum_mid(i) = vsum_mid(i) + (1._r8/gravit)*v(i,k)/rpdel(i,k)
              end do
           end do              

         
         
         
         
         
         
         
         
         
         
         

       else

         
         
         
         

         

           ksrf(:ncol) = ksrftms(:ncol) 

         

           u(:ncol,pver) = u(:ncol,pver) + tmp1(:ncol)*taux(:ncol)
           v(:ncol,pver) = v(:ncol,pver) + tmp1(:ncol)*tauy(:ncol)

       end if  

       
       
       
       
       
       
       
       
       

       call vd_lu_decomp( pcols , pver , ncol  ,                        &
                          ksrf  , kvm  , tmpi2 , rpdel , ztodt , zero , &
                          ca    , cc   , dnom  , tmpm  , ntop  , nbot )

       call vd_lu_solve(  pcols , pver , ncol  ,                        &
                          u     , ca   , tmpm  , dnom  , ntop  , nbot , zero )

       call vd_lu_solve(  pcols , pver , ncol  ,                        &
                          v     , ca   , tmpm  , dnom  , ntop  , nbot , zero )

       
       
       
       
       

       do i = 1, ncol

          
          
          
          

          tautmsx(i) = -ksrftms(i)*u(i,pver)
          tautmsy(i) = -ksrftms(i)*v(i,pver)

          if( do_iss ) then

            

              usum_out(i) = 0._r8
              vsum_out(i) = 0._r8
              do k = 1, pver
                 usum_out(i) = usum_out(i) + (1._r8/gravit)*u(i,k)/rpdel(i,k)
                 vsum_out(i) = vsum_out(i) + (1._r8/gravit)*v(i,k)/rpdel(i,k)
              end do

            
            
            
            
            

              tauimpx(i) = ( usum_out(i) - usum_in(i) ) / ztodt
              tauimpy(i) = ( vsum_out(i) - vsum_in(i) ) / ztodt

              tautotx(i) = tauimpx(i) 
              tautoty(i) = tauimpy(i) 

            
            
            

              if( itaures .eq. 1 ) then
                  tauresx(i) = taux(i) + tautmsx(i) + tauresx(i) - tauimpx(i)
                  tauresy(i) = tauy(i) + tautmsy(i) + tauresy(i) - tauimpy(i)
              endif

          else

              tautotx(i) = tautmsx(i) + taux(i)
              tautoty(i) = tautmsy(i) + tauy(i)
              tauresx(i) = 0._r8
              tauresy(i) = 0._r8

          end if  

       end do 

     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     

       
       
       

     
     

       
       
       
       

     
     

       k = pver + 1
       do i = 1, ncol
          tmpi1(i,1) = 0._r8
          tmpi1(i,k) = 0.5_r8 * ztodt * gravit * &
                       ( (-u(i,k-1) + dinp_u(i,k))*tautotx(i) + (-v(i,k-1) + dinp_v(i,k))*tautoty(i) )
       end do

       do k = 2, pver
          do i = 1, ncol
             dout_u = u(i,k) - u(i,k-1)
             dout_v = v(i,k) - v(i,k-1)
             tmpi1(i,k) = 0.25_r8 * tmpi2(i,k) * kvm(i,k) * &
                          ( dout_u**2 + dout_v**2 + dout_u*dinp_u(i,k) + dout_v*dinp_v(i,k) )
          end do
       end do

       

       do k = 1, pver
          do i = 1, ncol
             dtk(i,k) = ( tmpi1(i,k+1) + tmpi1(i,k) ) * rpdel(i,k)
             dse(i,k) = dse(i,k) + dtk(i,k)
          end do
       end do

    end if 

    
    
    

  
  

    if( diffuse(fieldlist,'s') ) then

      

        do k = 1, pver
           dse(:ncol,k) = dse(:ncol,k) + ztodt * rpdel(:ncol,k) * gravit  *                &
                                       ( rhoi(:ncol,k+1) * kvh(:ncol,k+1) * cgh(:ncol,k+1) &
                                       - rhoi(:ncol,k  ) * kvh(:ncol,k  ) * cgh(:ncol,k  ) )
       end do

     

       dse(:ncol,pver) = dse(:ncol,pver) + tmp1(:ncol) * shflx(:ncol)

     

       call vd_lu_decomp( pcols , pver , ncol  ,                         &
                          zero  , kvh  , tmpi2 , rpdel , ztodt , cc_top, &
                          ca    , cc   , dnom  , tmpm  , ntop  , nbot    )

       call vd_lu_solve(  pcols , pver , ncol  ,                         &
                          dse   , ca   , tmpm  , dnom  , ntop  , nbot  , cd_top )

     
     
     

       if( do_molec_diff ) then
           topflx(:ncol) =  - kvh(:ncol,ntop_molec) * tmpi2(:ncol,ntop_molec) / (ztodt*gravit) * &
                            ( dse(:ncol,ntop_molec) - dse_top(:ncol) )
       end if

    endif

    
    
    

  
  

    

    need_decomp = .true.

    do m = 1, ncnst

       if( diffuse(fieldlist,'q',m) ) then

           
           
           
           
           

           qtm(:ncol,:pver) = q(:ncol,:pver,m)

           do k = 1, pver
              q(:ncol,k,m) = q(:ncol,k,m) + &
                             ztodt * rpdel(:ncol,k) * gravit  * ( cflx(:ncol,m) * rrho(:ncol) ) * &
                           ( rhoi(:ncol,k+1) * kvh(:ncol,k+1) * cgs(:ncol,k+1)                    &
                           - rhoi(:ncol,k  ) * kvh(:ncol,k  ) * cgs(:ncol,k  ) )
           end do
           lqtst(:ncol) = all(q(:ncol,1:pver,m) >= qmincg(m), 2)
           do k = 1, pver
              q(:ncol,k,m) = merge( q(:ncol,k,m), qtm(:ncol,k), lqtst(:ncol) )
           end do

           

           q(:ncol,pver,m) = q(:ncol,pver,m) + tmp1(:ncol) * cflx(:ncol,m)

           

           if( need_decomp ) then

               call vd_lu_decomp( pcols , pver , ncol  ,                         &
                                  zero  , kvq  , tmpi2 , rpdel , ztodt , zero  , &
                                  ca    , cc   , dnom  , tmpm  , ntop  , nbot )

               if( do_molec_diff ) then

                 

                   status = vd_lu_qdecomp( pcols , pver   , ncol      , cnst_fixed_ubc(m), cnst_mw(m), ubc_mmr(:,m), &
                                           kvq   , kq_scal, mw_fac(m) , tmpi2            , rpdel     ,               &
                                           ca    , cc     , dnom      , tmpm             , rhoi      ,               &
                                           tint  , ztodt  , ntop_molec, nbot_molec       , cd_top )
               else
                   need_decomp =  .false.
               endif
           end if

           call vd_lu_solve(  pcols , pver , ncol  ,                         &
                              q(1,1,m) , ca, tmpm  , dnom  , ntop  , nbot  , cd_top )
       end if
    end do

    return
  end subroutine compute_vdiff

  
  
  

  subroutine vd_lu_decomp( pcols, pver, ncol ,                        &
                           ksrf , kv  , tmpi , rpdel, ztodt , cc_top, &
                           ca   , cc  , dnom , ze   , ntop  , nbot    )
    
    
    
    
    
    

    
    
    

    integer,  intent(in)  :: pcols                 
    integer,  intent(in)  :: pver                  
    integer,  intent(in)  :: ncol                  
    integer,  intent(in)  :: ntop                  
    integer,  intent(in)  :: nbot                  
    real(r8), intent(in)  :: ksrf(pcols)           
    real(r8), intent(in)  :: kv(pcols,pver+1)      
    real(r8), intent(in)  :: tmpi(pcols,pver+1)    
    real(r8), intent(in)  :: rpdel(pcols,pver)     
    real(r8), intent(in)  :: ztodt                 
    real(r8), intent(in)  :: cc_top(pcols)         

    real(r8), intent(out) :: ca(pcols,pver)        
    real(r8), intent(out) :: cc(pcols,pver)        
    real(r8), intent(out) :: dnom(pcols,pver)      
    real(r8), intent(out) :: ze(pcols,pver)        

    
    
    

    integer :: i                                   
    integer :: k                                   

    
    
    

    
    
    

    do k = nbot - 1, ntop, -1
       do i = 1, ncol
          ca(i,k  ) = kv(i,k+1) * tmpi(i,k+1) * rpdel(i,k  )
          cc(i,k+1) = kv(i,k+1) * tmpi(i,k+1) * rpdel(i,k+1)
       end do
    end do

    
    

    do i = 1, ncol
       ca(i,nbot) = 0._r8
    end do

    
    

    do i = 1, ncol
       dnom(i,nbot) = 1._r8/(1._r8 + cc(i,nbot) + ksrf(i)*ztodt*gravit*rpdel(i,nbot))
       ze(i,nbot)   = cc(i,nbot)*dnom(i,nbot)
    end do

    do k = nbot - 1, ntop + 1, -1
       do i = 1, ncol
          dnom(i,k) = 1._r8/(1._r8 + ca(i,k) + cc(i,k) - ca(i,k)*ze(i,k+1))
          ze(i,k)   = cc(i,k)*dnom(i,k)
       end do
    end do

    do i = 1, ncol
       dnom(i,ntop) = 1._r8/(1._r8 + ca(i,ntop) + cc_top(i) - ca(i,ntop)*ze(i,ntop+1))
    end do

    return
  end subroutine vd_lu_decomp

  
  
  

  subroutine vd_lu_solve( pcols , pver , ncol , &
                          q     , ca   , ze   , dnom , ntop , nbot , cd_top )
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    

    
    
    

    integer,  intent(in)    :: pcols                  
    integer,  intent(in)    :: pver                   
    integer,  intent(in)    :: ncol                   
    integer,  intent(in)    :: ntop                   
    integer,  intent(in)    :: nbot                   
    real(r8), intent(in)    :: ca(pcols,pver)         
    real(r8), intent(in)    :: ze(pcols,pver)         
    real(r8), intent(in)    :: dnom(pcols,pver)       
    real(r8), intent(in)    :: cd_top(pcols)          

    real(r8), intent(inout) :: q(pcols,pver)          

    
    
    

    real(r8)                :: zf(pcols,pver)         
    integer                    i, k                   

    
    
    

    
    
    

    do i = 1, ncol
       zf(i,nbot) = q(i,nbot)*dnom(i,nbot)
    end do

    do k = nbot - 1, ntop + 1, -1
       do i = 1, ncol
          zf(i,k) = (q(i,k) + ca(i,k)*zf(i,k+1))*dnom(i,k)
       end do
    end do

    

    k = ntop
    do i = 1, ncol
       zf(i,k) = (q(i,k) + cd_top(i) + ca(i,k)*zf(i,k+1))*dnom(i,k)
    end do

    

    do i = 1, ncol
       q(i,ntop) = zf(i,ntop)
    end do

    do k = ntop + 1, nbot, +1
       do i = 1, ncol
          q(i,k) = zf(i,k) + ze(i,k)*q(i,k-1)
       end do
    end do

    return
  end subroutine vd_lu_solve

  
  
  
  
  character(128) function vdiff_select( fieldlist, name, qindex )
    
    
    
    type(vdiff_selector), intent(inout)        :: fieldlist
    character(*),         intent(in)           :: name
    integer,              intent(in), optional :: qindex
    
    vdiff_select = ''
    select case (name)
    case ('u','U')
       fieldlist%fields(1) = .true.
    case ('v','V')
       fieldlist%fields(2) = .true.
    case ('s','S')
       fieldlist%fields(3) = .true.
    case ('q','Q')
       if( present(qindex) ) then
           fieldlist%fields(3 + qindex) = .true.
       else
           fieldlist%fields(4) = .true.
       endif
    case default
       write(vdiff_select,*) 'Bad argument to vdiff_index: ', name
    end select
    return
    
  end function vdiff_select

  type(vdiff_selector) function not(a)
    
    
    
    type(vdiff_selector), intent(in)  :: a
    allocate(not%fields(size(a%fields)))
    not%fields(:) = .not. a%fields(:)
  end function not

  logical function my_any(a)
    
    
    
    
    type(vdiff_selector), intent(in) :: a
    my_any = any(a%fields)
  end function my_any

  logical function diffuse(fieldlist,name,qindex)
    
    
    
    type(vdiff_selector), intent(in)           :: fieldlist
    character(*),         intent(in)           :: name
    integer,              intent(in), optional :: qindex
    
    select case (name)
    case ('u','U')
       diffuse = fieldlist%fields(1)
    case ('v','V')
       diffuse = fieldlist%fields(2)
    case ('s','S')
       diffuse = fieldlist%fields(3)
    case ('q','Q')
       if( present(qindex) ) then
           diffuse = fieldlist%fields(3 + qindex)
       else
           diffuse = fieldlist%fields(4)
       endif
    case default
       diffuse = .false.
    end select
    return
  end function diffuse

end module diffusion_solver
