module dmix ! ! (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland, ! P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED. ! implicit none complex, dimension(:,:), allocatable :: omegaz, omega0 real, dimension(:,:), allocatable :: e1, g1, dml, e_crit, D_old real, dimension(:,:), allocatable :: dm_out real, dimension(:), allocatable :: k_tim contains subroutine tc_out(g_target, d_target, err, icontrol, save_phi) ! ! Write critical temperature gradient info ! or D_R info ! use itg_data, only: dt, icrit, tol, ncycle, nstp, gamma_0, & etak, time, rdiff, runname, epsn, lrunname use fields, only: potential use mp, only: proc0, max_allreduce, barrier use gryffin_grid, only: ldb, r0, md, rky, rkperp2 use constants, only: ii use gryffin_layouts complex, dimension(:,m_low:,n_low:) :: save_phi real, dimension(:,:), allocatable :: e2, g2 complex w0,w1 complex phi1,phi0 real err real g_rate real D_max,g_max real g_target,d_target real dmix integer :: done = 0 integer i,l,m,n,init,icontrol(2),lsave,m_min if(.not.allocated(omegaz)) then allocate (omegaz(m_low:m_alloc, n_low:n_alloc), & e1(m_low:m_alloc, n_low:n_alloc), & D_old(m_low:m_alloc, n_low:n_alloc), & g1(m_low:m_alloc, n_low:n_alloc), & dml(m_low:m_alloc, n_low:n_alloc), & e_crit(m_low:m_alloc, n_low:n_alloc), & k_tim(m_low:m_alloc), & omega0(m_low:m_alloc, n_low:n_alloc), & dm_out(m_low:m_alloc, n_low:n_alloc)) endif m_min=2 if(md == 1) m_min=1 if(icontrol(2) == 2) init=1 if(done == 3) done=0 icontrol(2)=3 lsave=ldb/2+2 if(init == 1) then init=0 do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle D_old(m,n)=0. omegaz(m,n)=1. dml(m,n)=1. enddo enddo else g_max=-1.e10 do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle phi0=save_phi(lsave,m,n) phi1=potential(lsave,m,n) omega0(m,n)=omegaz(m,n) omegaz(m,n)=ii*log(phi1/phi0)/dt(m) g_max=max(g_max,aimag(omegaz(m,n))) D_old(m,n)=dml(m,n) enddo enddo call dm_calc(potential,omegaz,dml,D_max) endif if(icrit == 1) then err=0. do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle g_rate=aimag(omegaz(m,n)) w0=omega0(m,n) w1=omegaz(m,n) err=max(err,sqrt(cabs((w1-w0)/(w1+w0)))) err=max(err,abs((g_rate-g_target)/(g_target+g_rate))) enddo enddo call max_allreduce(err) endif if(icrit == 2) then err=0. do n=n_low, n_high do m=m_low, m_high if(m < m_min) cycle dmix=dml(m,n) err=max(err,abs((dmix-d_target)/(d_target+dmix))) err=max(err,abs((dmix-D_old(m,n))/(dmix+D_old(m,n)))) enddo enddo call max_allreduce(err) endif if(err < tol .or. ncycle == nstp) done=done+1 if(done == 1) then done=2 if(icontrol(1) == 0 .and. proc0) open(unit=1,file=runname(1:lrunname)//'.dmix') ! if(icrit == 1) write(*,*) D_max,g_max,etak(m_min,1,1), & ! ncycle,err,g_target,d_target if(icrit == 2) then ! write(*,*) D_max,ncycle,err,d_target ! write(*,*) (rdiff(j),j=1,n_ktheta) done=3 endif g_target=gamma_0/2. do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle g1(m,n)=aimag(omegaz(m,n)) e1(m,n)=etak(m,n,1) enddo enddo if(ncycle == nstp) done=3 endif if(done == 3) then if(icrit == 1) then ! write(*,*) D_max,g_max,etak(m_min,1,1),ncycle,err,g_target allocate (e2(m_low:m_alloc, n_low:n_alloc), g2(m_low:m_alloc, n_low:n_alloc)) do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle g2(m,n)=aimag(omegaz(m,n)) e2(m,n)=etak(m,n,1) if (g2(m,n)-g1(m,n) /= 0) then e_crit(m,n)=e2(m,n)-g2(m,n)*(e2(m,n)-e1(m,n))/(g2(m,n)-g1(m,n)) else e_crit(m,n)=-1.e10 end if enddo enddo deallocate (e2, g2) if(proc0) write(*,*) call barrier do n=n_low, n_high do m = m_low, m_high if(m < m_min) cycle if(e_crit(m,n) == -1.e10) then write(*,2002) rky(m,n),r0(m,n) else write(*,2001) time(m),rky(m,n),r0(m,n),& e_crit(m,n),e_crit(m,n)/epsn endif enddo enddo call barrier if(proc0) write(*,*) call barrier icontrol(1)=1 endif if(icrit == 2) then do n = n_low, n_high do m = m_low, m_high if (m < m_min) cycle dm_out(m,n)=dml(m,n)+rdiff(m,n) enddo enddo call barrier if(proc0) write(*,*) call barrier do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle write(*,2000) time(m),rky(m,n),r0(m,n),dm_out(m,n) enddo enddo call barrier if(proc0) write(*,*) call barrier icontrol(1)=1 endif return endif 2000 format(' time= ',f7.2,' k_theta= ',f6.4,' theta_0= ',f8.6,' D_ML= ',f8.5) 2001 format(' time= ',f7.2,' k_theta= ',f6.4,' theta_0= ',f6.4,' Ln/LT_crit= ',f8.5,' R/LT_crit= ',f8.5) 2002 format('Failed to find LT_crit for k_theta= ',f6.4,' theta_0= ',f6.4) end subroutine tc_out subroutine dm_calc(a, omegaz, dml, D_max) use gryffin_grid, only: ldb, rkperp2, md use gryffin_layouts complex, dimension(:,m_low:,n_low:) :: a complex, dimension(m_low:,n_low:) :: omegaz real, dimension(m_low:,n_low:) :: dml real :: norm, D_max, avekp2 integer :: l, m, n D_max=0. do n = n_low, n_high do m = m_low, m_high avekp2=0.0 norm=0.0 do l=1,ldb avekp2=avekp2+rkperp2(l,m,n)*cabs(a(l,m,n))**2 norm=norm+cabs(a(l,m,n))**2 enddo if (norm /= 0.0 .and. avekp2 /= 0.) then avekp2=avekp2/norm dml(m,n)=aimag(omegaz(m,n))/avekp2 else dml(m,n)=0. endif D_max=max(D_max,dml(m,n)) enddo enddo end subroutine dm_calc subroutine critical(g_target,d_target,icontrol2) ! ! Find the critical temperature gradient or the diffusion needed ! to stabilize the mode. ! ! gamma_0 should be made a function of k_theta. ! use itg_data, only: etak, etak_par, rdiff, rdiff_0, time, & icrit, nspecies, gamma_0 use fields, only: potential use gryffin_grid, only: ldb, rkperp2, md use gryffin_layouts use diagnostics, only: volsq use linear, only: flrinit real, dimension(:,:), pointer, save :: p2, etak_0, d_0, k_p2 real g_target, d_target real dum, d_fac, e_fac, avekp2 integer i, l, m, n, init, icontrol2, m_min save init if(.not.associated(p2)) then allocate (p2(m_low:m_alloc, n_low:n_alloc), & etak_0(m_low:m_alloc, n_low:n_alloc), & d_0(m_low:m_alloc, n_low:n_alloc), & k_p2(m_low:m_alloc, n_low:n_alloc)) endif ! ! Assume m=1 is ky=0 mode unless md=1; don't calculate anything for ! ky=0 mode here: ! m_min=2 if(md == 1) m_min=1 if(icontrol2 == 1) then init=icontrol2 icontrol2=2 endif if(init == 1) then init=2 g_target=gamma_0 d_target=rdiff_0 call volsq(potential,potential,dum,p2) do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle etak_0(m,n)=etak(m,n,1) d_0(m,n)=rdiff(m,n) k_tim(m)=time(m) k_p2(m,n)=p2(m,n) enddo enddo else if(init == 2) then call volsq(potential,potential,dum,p2) if(icrit == 1) then if(gamma_0 /= g_target) then do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle etak_0(m,n)=etak(m,n,1) k_tim(m)=time(m) k_p2(m,n)=p2(m,n) enddo enddo gamma_0=g_target endif do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle e_fac=exp(2.*gamma_0*(time(m)-k_tim(m))) ! write(*,*) e_fac, etak_0(m,:) ! write(*,*) ! write(*,*) p2(m,:) ! write(*,*) ! write(*,*) k_p2(m,:) ! write(*,*) 'end' etak(m,n,1)=etak_0(m,n)/sqrt(p2(m,n)/k_p2(m,n)/e_fac) etak_par(m,n,1)=etak(m,n,1) enddo enddo if(nspecies > 1) then do i=2,nspecies do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle etak(m,n,i)=etak(m,n,1) etak_par(m,n,i)=etak_par(m,n,1) enddo enddo enddo endif call flrinit(0) else if(icrit == 2) then if(rdiff_0 /= d_target) then do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle d_0(m,n)=rdiff(m,n) k_tim(m)=time(m) k_p2(m,n)=p2(m,n) enddo enddo rdiff_0=d_target endif do n = n_low, n_high do m = m_low, m_high if(m < m_min) cycle avekp2=0. dum=0. do l=1,ldb avekp2=avekp2+cabs(potential(l,m,n))**2*rkperp2(l,m,n) dum=dum+cabs(potential(l,m,n))**2 enddo if(dum /= 0.) then avekp2=avekp2/dum else avekp2=0. endif d_fac=exp(2.*rdiff_0*avekp2*(time(m)-k_tim(m))) rdiff(m,n)=d_0(m,n)*sqrt(p2(m,n)/k_p2(m,n)/d_fac) enddo enddo endif else return endif end subroutine critical end module dmix