EQUIL - FEATURE - Optional profile_source = "integrate": build F and P in the sq spline from FF′ and p′ (values stays the default)

docs/examples/DIIID_resistive_example/rdcon.in

@@ -75,8 +75,8 @@
     bin_bal1=f           ! Binary output for bal_flag poloidal functions
     out_bal2=f           ! Ascii output for bal_flag functions
     bin_bal2=f           ! Binary output for bal_flag functions
-    MRE_flag=t           ! Netcdf output for various modified Rutherford equation terms (defined in mercier.f)
-    geom_flag=t          ! Netcdf output for more flux-surface averaged terms (defined in mercier.f). To use, must set MRE_flag = t
+    MRE_flag=f           ! Netcdf output for various modified Rutherford equation terms (defined in mercier.f)
+    geom_flag=f          ! Netcdf output for more flux-surface averaged terms (defined in mercier.f). To use, must set MRE_flag = t
 /
 
 &UA_DIAGNOSE_LIST

equil/direct.f

@@ -74,7 +74,7 @@ SUBROUTINE direct_run
       sq_in%fs(:,4)=SQRT(sq_in%xs)
       sq_in%name="  sq  "
       sq_in%title=(/"psifac","  f   ","mu0 p ","  q   "," rho  "/)
-      CALL spline_fit(sq_in,"extrap")
+      CALL spline_fit_hermite(sq_in,"extrap",sq_in_slopes)
       psi_in%xs=rmin+(/(ir,ir=0,psi_in%mx)/)*(rmax-rmin)/psi_in%mx
       psi_in%ys=zmin+(/(iz,iz=0,psi_in%my)/)*(zmax-zmin)/psi_in%my
       CALL bicube_fit(psi_in,"extrap","extrap")
@@ -179,6 +179,8 @@ SUBROUTINE direct_run
 c-----------------------------------------------------------------------
          sq%fs(ipsi,1)=bf%f*twopi
          sq%fs(ipsi,2)=bf%p
+         sq%fs1(ipsi,1)=bf%f1*twopi
+         sq%fs1(ipsi,2)=bf%p1
          sq%fs(ipsi,3)=y_out(istep,1)*twopi*psio
          sq%fs(ipsi,4)=y_out(istep,3)*bf%f/twopi
          CALL spline_dealloc(ff)
@@ -191,9 +193,8 @@ SUBROUTINE direct_run
 c-----------------------------------------------------------------------
 c     fit surface quantities to cubic splines.
 c-----------------------------------------------------------------------
-      CALL spline_fit(sq,"extrap")
+      CALL spline_fit_hermite(sq,"extrap",sq_in_slopes)
       sq%name="  sq  "
-      sq%title=(/" psi  ","  f   ","  p   ","  q   "/)
       q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0)
 c-----------------------------------------------------------------------
 c     revise q profile.
@@ -205,10 +206,11 @@ SUBROUTINE direct_run
          DO ipsi=0,mpsi
             ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)
             sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac
+            sq%fs1(ipsi,1)=sq%fs1(ipsi,1)/ffac
             sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac
             rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac
          ENDDO
-         CALL spline_fit(sq,"extrap")
+         CALL spline_fit_hermite(sq,"extrap",sq_in_slopes)
       ENDIF
       qa=sq%fs(mpsi,4)+sq%fs1(mpsi,4)*(1-sq%xs(mpsi))
 c-----------------------------------------------------------------------

equil/equil.f

@@ -39,7 +39,7 @@ SUBROUTINE equil_loadnamelists(op_psihigh, op_psilow)
      $     newq0,psihigh,psilow,input_only,jac_type,power_bp,power_r,
      $     power_b,jac_method,convert_type,power_flag,ns1,
      $     sp_pfac,sp_nx,sp_dx1,sp_dx2,use_galgrid,etol,enstep,
-     $     use_classic_splines
+     $     use_classic_splines,profile_source
       NAMELIST/equil_output/bin_2d,bin_eq_1d,bin_eq_2d,out_2d,out_eq_1d,
      $     out_eq_2d,bin_fl,out_fl,interp,gse_flag,dump_flag,verbose,
      $     out_ahg2msc
@@ -150,6 +150,7 @@ SUBROUTINE equil_read(unit, op_psihigh, op_psilow)
 c-----------------------------------------------------------------------
 c     read equilibrium data and diagnose.
 c-----------------------------------------------------------------------
+      sq_in_slopes=.FALSE.
       SELECT CASE(TRIM(eq_type))
       CASE("fluxgrid")
          CALL read_eq_fluxgrid

equil/equil_out.f

@@ -97,7 +97,9 @@ SUBROUTINE equil_out_diagnose(sq_flag,unit)
          sq_out%fs(:,2)=sq%fs(:,2)
          sq_out%fs(:,3)=sq%fs(:,4)
          sq_out%fs(:,4)=SQRT(sq_out%xs)
-         CALL spline_fit(sq_out,"extrap")
+         sq_out%fs1(:,1)=sq%fs1(:,1)/twopi
+         sq_out%fs1(:,2)=sq%fs1(:,2)
+         CALL spline_fit_hermite(sq_out,"extrap",sq_in_slopes)
          sq_out%name="  sq  "
          sq_out%title=(/"psifac","  f   ","mu0 p ","  q   "," rho  "/)
          IF(bin_eq_1d)CALL bin_open(bin_2d_unit,"sq_out.bin","UNKNOWN",
@@ -845,6 +847,11 @@ SUBROUTINE equil_out_dump
       WRITE(dump_unit)rzphi%name
       WRITE(dump_unit)rzphi%periodic
 c-----------------------------------------------------------------------
+c     write f, p knot derivatives and eqfun; optional for read_eq_dump.
+c-----------------------------------------------------------------------
+      WRITE(dump_unit)sq_in_slopes(1:2),sq%fs1(:,1:2)
+      IF(eqfun%allocated)WRITE(dump_unit)eqfun%fs
+c-----------------------------------------------------------------------
 c     close dump file.
 c-----------------------------------------------------------------------
       CALL bin_close(dump_unit)

equil/global.f

@@ -16,6 +16,8 @@ MODULE global_mod
 
       CHARACTER(16) :: eq_type="fluxgrid"
       CHARACTER(128) :: eq_filename="fluxgrid.dat"
+      CHARACTER(16) :: profile_source="values"
+      LOGICAL, DIMENSION(4) :: sq_in_slopes=.FALSE.
 
       LOGICAL :: direct_flag=.FALSE.
       LOGICAL :: out_eq_1d=.FALSE.,bin_eq_1d=.FALSE.

equil/inverse.f

@@ -54,7 +54,7 @@ SUBROUTINE inverse_run
       sq_in%fs(:,4)=SQRT(sq_in%xs)
       sq_in%name="  sq  "
       sq_in%title=(/"psifac","  f   ","mu0 p ","  q   "," rho  "/)
-      CALL spline_fit(sq_in,"extrap")
+      CALL spline_fit_hermite(sq_in,"extrap",sq_in_slopes)
       rz_in%xs=sq_in%xs
       rz_in%ys=(/(itheta,itheta=0,rz_in%my)/)/REAL(rz_in%my,r8)
       rz_in%name="  rz  "
@@ -180,7 +180,7 @@ SUBROUTINE inverse_run
       CALL spline_alloc(spl,mtheta,5)
       DO ipsi=0,mpsi
          psifac=rzphi%xs(ipsi)
-         CALL spline_eval(sq_in,psifac,0)
+         CALL spline_eval(sq_in,psifac,1)
          spl%xs=rzphi%ys
 c-----------------------------------------------------------------------
 c     interpolate to new surface.
@@ -232,10 +232,12 @@ SUBROUTINE inverse_run
 c-----------------------------------------------------------------------
          sq%fs(ipsi,1)=sq_in%f(1)*twopi
          sq%fs(ipsi,2)=sq_in%f(2)
+         sq%fs1(ipsi,1)=sq_in%f1(1)*twopi
+         sq%fs1(ipsi,2)=sq_in%f1(2)
          sq%fs(ipsi,3)=spl%fsi(mtheta,3)*twopi*pi
          sq%fs(ipsi,4)=spl%fsi(mtheta,4)*sq%fs(ipsi,1)/(2*twopi*psio)
       ENDDO
-      CALL spline_fit(sq,"extrap")
+      CALL spline_fit_hermite(sq,"extrap",sq_in_slopes)
       q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0)
 c-----------------------------------------------------------------------
 c     revise q profile.
@@ -247,10 +249,11 @@ SUBROUTINE inverse_run
          DO ipsi=0,mpsi
             ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)
             sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac
+            sq%fs1(ipsi,1)=sq%fs1(ipsi,1)/ffac
             sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac
             rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac
          ENDDO
-         CALL spline_fit(sq,"extrap")
+         CALL spline_fit_hermite(sq,"extrap",sq_in_slopes)
       ENDIF
       qa=sq%fs(mpsi,4)+sq%fs1(mpsi,4)*(one-sq%xs(mpsi))
 c-----------------------------------------------------------------------

equil/local.f

@@ -13,6 +13,7 @@
 c     5. ascii_close.
 c     6. program_stop.
 c     7. floored_log.
+c     8. read_var_len.
 c-----------------------------------------------------------------------
 c     subprogram 0. local.
 c     module declarations.
@@ -259,4 +260,46 @@ FUNCTION floored_log(u) RESULT(ulog)
 c-----------------------------------------------------------------------
       RETURN
       END FUNCTION floored_log
+c-----------------------------------------------------------------------
+c     subprogram 8. read_var_len.
+c     allocates variable length arrays from set-length read-in arrays.
+c-----------------------------------------------------------------------
+c-----------------------------------------------------------------------
+c     declarations.
+c-----------------------------------------------------------------------
+      SUBROUTINE read_var_len(x_in, x_out)
+
+      IMPLICIT NONE
+
+      REAL(r8), DIMENSION(:), INTENT(IN) :: x_in ! Indexes from 1...
+      REAL(r8), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: x_out
+      INTEGER :: n, i
+      REAL(r8) :: default_value
+c-----------------------------------------------------------------------
+c     assume last value of x_in is a sentinal
+c-----------------------------------------------------------------------
+      default_value = x_in(SIZE(x_in))
+c-----------------------------------------------------------------------
+c     count leading non-sentinel values in x_in.
+c-----------------------------------------------------------------------
+      n = 0
+      DO i = 1, SIZE(x_in) - 1
+        IF (x_in(i) /= default_value) THEN
+          n = n + 1
+        ELSE
+          EXIT
+        END IF
+      END DO
+c-----------------------------------------------------------------------
+c     allocate output array and copy if values exist.
+c-----------------------------------------------------------------------
+      IF (n > 0) THEN
+        ALLOCATE(x_out(n))
+        x_out(1:n) = x_in(1:n)
+      END IF
+c-----------------------------------------------------------------------
+c     terminate.
+c-----------------------------------------------------------------------
+      RETURN
+      END SUBROUTINE read_var_len
       END MODULE local_mod