Many updates for new global shearing terms

cgyro/src/cgyro_advect_wavenumber.F90

@@ -2,7 +2,7 @@
 ! cgyro_advect_wavenumber.f90
 !
 ! PURPOSE:
-!  Manage shearing by wavenumber advection.
+!  Manage shearing by wavenumber advection (legacy method)
 !---------------------------------------------------------
 
 subroutine cgyro_advect_wavenumber(ij)

cgyro/src/cgyro_globalshear.F90

@@ -1,9 +1,9 @@
-!---------------------------------------------------------
-! cgyro_advect_wavenumber.f90
+!-----------------------------------------------------------------
+! cgyro_globalshear.F90
 !
 ! PURPOSE:
-!  Manage shearing by wavenumber advection.
-!---------------------------------------------------------
+!  Manage shearing by wavenumber advection (new Dudkovskaia terms)
+!-----------------------------------------------------------------
 
 subroutine cgyro_globalshear(ij)
 
@@ -47,16 +47,22 @@ subroutine cgyro_globalshear(ij)
                  llnt = ll*n_theta
 
                  if ( (ir+ll) <= n_radial ) then
+                    ! ExB shear
                     h1 = omega_eb_base*itor*h_x(iccj+llnt,ivc,itor)
+                    ! omega_star shear
                     h1 = h1-sum(omega_ss(:,iccj+llnt,ivc,itor)*field(:,iccj+llnt,itor))
+                    ! beta_star shear
                     h1 = h1-omega_sbeta(iccj+llnt,ivc,itor)*cap_h_c(iccj+llnt,ivc,itor)
                  else
                     h1 = 0.0
                  endif
 
                  if ( (ir-ll) >= 1 ) then
+                    ! ExB shear
                     h2 = omega_eb_base*itor*h_x(iccj-llnt,ivc,itor)
+                    ! omega_star shear
                     h2 = h2-sum(omega_ss(:,iccj-llnt,ivc,itor)*field(:,iccj-llnt,itor))
+                    ! beta_star shear
                     h2 = h2-omega_sbeta(iccj-llnt,ivc,itor)*cap_h_c(iccj-llnt,ivc,itor)
                  else
                     h2 = 0.0

cgyro/src/cgyro_init_arrays.F90

@@ -456,15 +456,17 @@ subroutine cgyro_init_arrays
            ! global-Taylor corrections via wavenumber advection (ix -> d/dp)
            ! JC: Re-checked sign and normalization (Oct 2019)
 
-           ! generalized profile curvature (phi shearing)
+           ! NEW GLOBAL TERMS
+
+           ! generalized omega_star shear (acts on phi)
            sm = sdlnndr(is)+sdlntdr(is)*(energy(ie)-1.5)
 
-           ! beta_star (H shearing)
-           sb = sbeta_star(is)*energy(ie)/bmag(it)**2
+           ! generalized beta/drift shear (acts on H)
+           sb = beta_star(0)*sbeta_star(is)*energy(ie)/bmag(it)**2
 
            arg = -k_theta_base*itor*length/(2*pi)
 
-           omega_ss(:,ic,iv_loc,itor) = arg*sm*jvec_c(:,ic,iv_loc,itor)
+           omega_ss(:,ic,iv_loc,itor)  = arg*sm*jvec_c(:,ic,iv_loc,itor)
            omega_sbeta(ic,iv_loc,itor) = arg*sb
 
         enddo

f2py/expro/expro.f90

@@ -121,6 +121,7 @@ module expro
        expro_dlntedr,&
        expro_sdlnnedr,&
        expro_sdlntedr,&
+       expro_sbetae,&
        expro_dlnptotdr,&
        expro_w0p,&
        expro_surf,&
@@ -160,7 +161,8 @@ module expro
        expro_dlnnidr,&
        expro_dlntidr,&
        expro_sdlnnidr,&
-       expro_sdlntidr
+       expro_sdlntidr,&
+       expro_sbetai
   !-----------------------------------------------------------
 
   ! input.gacode.geo dimension and arrays
@@ -287,6 +289,7 @@ subroutine expro_init(flag)
        allocate(expro_dlntedr(nexp))      ; expro_dlntedr = 0.0
        allocate(expro_sdlnnedr(nexp))     ; expro_sdlnnedr = 0.0
        allocate(expro_sdlntedr(nexp))     ; expro_sdlntedr = 0.0
+       allocate(expro_sbetae(nexp))       ; expro_sbetae = 0.0
        allocate(expro_dlnptotdr(nexp))    ; expro_dlnptotdr = 0.0
        allocate(expro_w0p(nexp))          ; expro_w0p = 0.0
        allocate(expro_surf(nexp))         ; expro_surf = 0.0
@@ -328,6 +331,7 @@ subroutine expro_init(flag)
        allocate(expro_dlntidr(nion,nexp))  ; expro_dlntidr = 0.0
        allocate(expro_sdlnnidr(nion,nexp)) ; expro_sdlnnidr = 0.0
        allocate(expro_sdlntidr(nion,nexp)) ; expro_sdlntidr = 0.0
+       allocate(expro_sbetai(nion,nexp))   ; expro_sbetai = 0.0
 
     else
 
@@ -451,6 +455,7 @@ subroutine expro_init(flag)
        deallocate(expro_dlntidr)  
        deallocate(expro_sdlnnidr) 
        deallocate(expro_sdlntidr)  
+       deallocate(expro_sbetai)  
 
     endif
 

f2py/expro/expro_locsim.f90

@@ -14,6 +14,7 @@ module expro_locsim_interface
   double precision, dimension(:,:), allocatable :: dlnndr_exp
   double precision, dimension(:,:), allocatable :: sdlntdr_exp
   double precision, dimension(:,:), allocatable :: sdlnndr_exp
+  double precision, dimension(:,:), allocatable :: sbeta_exp
 
   double precision, dimension(:), allocatable :: gamma_e_exp
   double precision, dimension(:), allocatable :: gamma_p_exp
@@ -79,6 +80,7 @@ module expro_locsim_interface
   double precision, dimension(9) :: dlntdr_loc
   double precision, dimension(9) :: sdlnndr_loc
   double precision, dimension(9) :: sdlntdr_loc
+  double precision, dimension(9) :: sbeta_loc
 
   integer :: geo_ny_loc
   double precision, dimension(:,:), allocatable :: geo_yin_loc
@@ -106,6 +108,7 @@ subroutine expro_locsim_alloc(flag)
      allocate(dlnndr_exp(n_species_exp,expro_n_exp))
      allocate(sdlntdr_exp(n_species_exp,expro_n_exp))
      allocate(sdlnndr_exp(n_species_exp,expro_n_exp))
+     allocate(sbeta_exp(n_species_exp,expro_n_exp))
 
      allocate(gamma_e_exp(expro_n_exp))
      allocate(gamma_p_exp(expro_n_exp))
@@ -121,6 +124,7 @@ subroutine expro_locsim_alloc(flag)
      if (allocated(dlnndr_exp)) deallocate(dlnndr_exp)
      if (allocated(sdlntdr_exp)) deallocate(sdlntdr_exp)
      if (allocated(sdlnndr_exp)) deallocate(sdlnndr_exp)
+     if (allocated(sbeta_exp)) deallocate(sbeta_exp)
 
      if (allocated(gamma_e_exp)) deallocate(gamma_e_exp)
      if (allocated(gamma_p_exp)) deallocate(gamma_p_exp)
@@ -279,6 +283,7 @@ subroutine expro_locsim_profiles(&
   dens_exp(n_species_exp,:)    = expro_ne(:)
   dlnndr_exp(n_species_exp,:)  = expro_dlnnedr(:)*a_meters 
   sdlnndr_exp(n_species_exp,:) = expro_sdlnnedr(:)*a_meters
+  sbeta_exp(n_species_exp,:)   = expro_sbetae(:)*a_meters
 
   mass_loc(n_species_exp) = expro_masse
   z_loc(n_species_exp) = -1d0
@@ -303,6 +308,9 @@ subroutine expro_locsim_profiles(&
         dlnndr_exp(i_ion,:)  = expro_dlnnidr(i_ion,:)*a_meters
         sdlnndr_exp(i_ion,:) = expro_sdlnnidr(i_ion,:)*a_meters
      endif
+
+     sbeta_exp(i_ion,:) = expro_sbetai(i_ion,:)*a_meters 
+
   enddo
 
   ! Rotation
@@ -369,11 +377,12 @@ subroutine expro_locsim_profiles(&
      call cub_spline1(rmin_exp,dlnndr_exp(i,:),expro_n_exp,rmin,dlnndr_loc(i))
      call cub_spline1(rmin_exp,sdlntdr_exp(i,:),expro_n_exp,rmin,sdlntdr_loc(i))
      call cub_spline1(rmin_exp,sdlnndr_exp(i,:),expro_n_exp,rmin,sdlnndr_loc(i))
+     call cub_spline1(rmin_exp,sbeta_exp(i,:),expro_n_exp,rmin,sbeta_loc(i))
      beta_star_loc = beta_star_loc+dens_loc(i)*temp_loc(i)*(dlnndr_loc(i)+dlntdr_loc(i))
   enddo
-  ! CGS beta calculation
+
+  ! CGS beta calculation (JC: is this needed?)
   betae_loc = 4.027e-3*dens_loc(n_species_exp)*temp_loc(n_species_exp)/b_unit_loc**2
-
   beta_star_loc = beta_star_loc*betae_loc/(dens_loc(n_species_exp)*temp_loc(n_species_exp))
 
   if (numeq_flag == 1 .and. expro_nfourier > 0) then

f2py/expro/expro_util.f90

@@ -22,6 +22,8 @@ subroutine expro_compute_derived
   double precision, dimension(:), allocatable :: temp
   double precision, dimension(:), allocatable :: cc
   double precision, dimension(:), allocatable :: loglam
+  double precision, dimension(:,:), allocatable :: sni,sti
+  double precision, dimension(:), allocatable  :: sne,ste
 
   double precision :: r_min
   double precision :: fa,fb
@@ -30,6 +32,7 @@ subroutine expro_compute_derived
   double precision :: eps0,m0,ipma
   double precision :: bt2_ave
   double precision :: bp2_ave
+
 
   mp = expro_mass_deuterium/2d0  ! mass_deuterium/2.0 (g)
 
@@ -64,6 +67,10 @@ subroutine expro_compute_derived
   !
   allocate(torflux(expro_n_exp))
   allocate(temp(expro_n_exp))
+  allocate(sne(expro_n_exp))
+  allocate(ste(expro_n_exp))
+  allocate(sni(expro_n_ion,expro_n_exp))
+  allocate(sti(expro_n_ion,expro_n_exp))
 
   torflux(:) = expro_torfluxa*expro_rho(:)**2
 
@@ -136,19 +143,19 @@ subroutine expro_compute_derived
   call bound_deriv(expro_dlntedr,-log(expro_te),expro_rmin,expro_n_exp)
 
   ! NOTE: expro_sdln* will be renormalized after calculation of rhos later
-
+  
   ! sne = -ne''/ne (1/m^2) [not fully normalized yet]
-  call bound_deriv(expro_sdlnnedr,expro_ne*expro_dlnnedr,expro_rmin,expro_n_exp)
-  expro_sdlnnedr = expro_sdlnnedr/expro_ne
+  call bound_deriv(sne,expro_ne*expro_dlnnedr,expro_rmin,expro_n_exp)
+  sne = sne/expro_ne
 
   ! sTe = -Te''/Te (1/m^2) [not fully normalized yet]
-  call bound_deriv(expro_sdlntedr,expro_te*expro_dlntedr,expro_rmin,expro_n_exp)
-  expro_sdlntedr = expro_sdlntedr/expro_te
+  call bound_deriv(ste,expro_te*expro_dlntedr,expro_rmin,expro_n_exp)
+  ste = ste/expro_te
 
   expro_dlnnidr = 0d0
   expro_dlntidr = 0d0
-  expro_sdlnnidr = 0d0
-  expro_sdlntidr = 0d0
+  sni = 0d0
+  sti = 0d0
 
   do is=1,expro_n_ion
      if (minval(expro_ni(is,:)) > 0d0 .and. minval(expro_ti(is,:)) > 0d0) then
@@ -159,21 +166,37 @@ subroutine expro_compute_derived
         call bound_deriv(expro_dlntidr(is,:),-log(expro_ti(is,:)),expro_rmin,expro_n_exp)
 
         ! sni = -ni''/ni (1/m^2) [not fully normalized yet]
-        call bound_deriv(expro_sdlnnidr(is,:),expro_ni(is,:)*expro_dlnnidr(is,:),expro_rmin,expro_n_exp)
-        expro_sdlnnidr(is,:) = expro_sdlnnidr(is,:)/expro_ni(is,:)
+        call bound_deriv(sni(is,:),expro_ni(is,:)*expro_dlnnidr(is,:),expro_rmin,expro_n_exp)
+        sni(is,:) = sni(is,:)/expro_ni(is,:)
 
         ! sTi = -Ti''/Ti (1/m^2) [not fully normalized yet]
         call bound_deriv(expro_sdlntidr(is,:),expro_ti(is,:)*expro_dlntidr(is,:),expro_rmin,expro_n_exp)
-        expro_sdlntidr(is,:) = expro_sdlntidr(is,:)/expro_ti(is,:)
+        sti(is,:) = sti(is,:)/expro_ti(is,:)
      endif
   enddo
 
+  ! omega_star shear
+  expro_sdlnnedr = sne+1.5*expro_dlntedr**2+expro_dlnnedr**2-expro_dlnnedr*expro_dlntedr
+  expro_sdlnnidr = sni-1.5*expro_dlntidr**2+expro_dlnnidr**2-expro_dlnnidr*expro_dlntidr
+
+  expro_sdlntedr = ste+expro_dlntedr**2
+  expro_sdlntidr = sti+expro_dlntidr**2
+
+  ! beta_star shear 
+  expro_sbetae = (sne+ste-2*expro_dlntedr*expro_dlnnedr)/(expro_dlnnedr+expro_dlntedr)
+  expro_sbetai = (sni+sti-2*expro_dlntidr*expro_dlnnidr)/(expro_dlnnidr+expro_dlntidr)
+
   ! 1/L_Ptot = -dln(Ptot)/dr (1/m)
   if (minval(expro_ptot) > 0d0) then
      call bound_deriv(expro_dlnptotdr,-log(expro_ptot),expro_rmin,expro_n_exp)
   else
      expro_dlnptotdr = 0d0
   endif
+
+  deallocate(sne)
+  deallocate(ste)
+  deallocate(sni)
+  deallocate(sti)
   !--------------------------------------------------------------------
 
   !-------------------------------------------------------------------
@@ -314,7 +337,6 @@ subroutine expro_compute_derived
   expro_vol(1)  = 0d0
   expro_volp(1) = 0d0  
   expro_thetascale(1) = expro_thetascale(2)
-
   !--------------------------------------------------------------
 
   !-----------------------------------------------------------------
@@ -335,11 +357,14 @@ subroutine expro_compute_derived
   !
   ! sn = -n''/n*rhos (1/m) 
   ! sT = -T''/T*rhos (1/m)
+  ! sbeta = sbeta*rhos (1/m)
   expro_sdlnnedr = expro_sdlnnedr*expro_rhos(:)
   expro_sdlntedr = expro_sdlntedr*expro_rhos(:)
+  expro_sbetae = expro_sbetae*expro_rhos(:)
   do is=1,expro_n_ion
      expro_sdlnnidr(is,:) = expro_sdlnnidr(is,:)*expro_rhos(:)
      expro_sdlntidr(is,:) = expro_sdlntidr(is,:)*expro_rhos(:)
+     expro_sbetai(is,:) = expro_sbetai(is,:)*expro_rhos(:)
   enddo
   !-----------------------------------------------------------------
 

profiles_gen/locpargen/locpargen_cgyro.f90

@@ -113,6 +113,7 @@ subroutine locpargen_cgyro
      write(1,10) 'DLNTDR_'//tag(is)//'=',dlntdr_loc(is)
      write(1,10) 'SDLNNDR_'//tag(is)//'=',sdlnndr_loc(is)
      write(1,10) 'SDLNTDR_'//tag(is)//'=',sdlntdr_loc(is)
+     write(1,10) 'SBETA_'//tag(is)//'=',sbeta_loc(is)
   enddo
   close(1)