anchor.diff:               do k = 1, ADM_kall
                do g = 1, ADM_gall_pl
                   d2T_dx2 = abs(vtmp_pl(g,k,l,5)) / T0 * AREA_ave
  -                coef    = cfact * ( AREA_ave * AREA_ave ) / TIME_DTL * d2T_dx2
  +                coef    = cfact * ( AREA_ave * AREA_ave ) / large_step_dt * d2T_dx2
   
                   KH_coef_pl(g,k,l) = max( min( coef, Kh_max(k) ), Kh_coef_minlim )
                enddo
  
positive.code:        if ( p == lap_order_hdiff ) then
  
            if ( hdiff_nonlinear ) then
               do l = 1, ADM_lall
               do k = 1, ADM_kall
               do g = 1, ADM_gall
                  d2T_dx2 = abs(vtmp(g,k,l,5)) / T0 * AREA_ave
                  coef    = cfact * ( AREA_ave * AREA_ave ) / TIME_DTL * d2T_dx2
  
                  KH_coef(g,k,l) = max( min( coef, Kh_max(k) ), Kh_coef_minlim )
               enddo
  
negative.code:     use mod_gtl, only: &
         GTL_max_k, &
         GTL_min_k
      use mod_runconf, only: &
         DYN_DIV_NUM
      implicit none
  
      character(len=*), intent(in) :: divdamp_type ! type of divergence damping
  