A new branch for snu students to incorporate local developments into …

dcon/dcon.F

@@ -22,6 +22,7 @@ PROGRAM dcon
       USE ode_mod
       USE free_mod
       USE resist_mod
+      USE bernstein_mod
       USE pentrc_interface,         ! rename overlapping names
      $    pentrc_verbose=>verbose,  ! should get a more fundamental fix
      $    pentrc_mpert=>mpert,
@@ -66,7 +67,7 @@ PROGRAM dcon
      $     out_sol_min,out_sol_max,bin_sol,bin_sol_min,bin_sol_max,
      $     out_fl,bin_fl,out_evals,bin_evals,bin_euler,euler_stride,
      $     bin_vac,ahb_flag,mthsurf0,msol_ahb,netcdf_out,out_fund,
-     $     out_ahg2msc
+     $     out_ahg2msc,bernstien_flag
 c-----------------------------------------------------------------------
 c     format statements.
 c-----------------------------------------------------------------------
@@ -136,6 +137,9 @@ PROGRAM dcon
 c-----------------------------------------------------------------------
       CALL equil_out_diagnose(.FALSE.,out_unit)
       CALL equil_out_write_2d
+      IF (bernstien_flag) THEN
+          CALL bernstein_calculation
+      END IF
       IF(dump_flag .AND. eq_type /= "dump")CALL equil_out_dump
       IF(direct_flag)CALL bicube_dealloc(psi_in)
 c-----------------------------------------------------------------------
@@ -238,7 +242,7 @@ PROGRAM dcon
             CALL read_kin(kinetic_file,zi,zimp,mi,mimp,nfac,
      $          tfac,wefac,wpfac,indebug)
             ! manually set the perturbed equilibrium displacements
-            ! use false flat xi and xi' for equal weighting
+            ! use false flat xi and xi for equal weighting
             ALLOCATE(psitmp(sq%mx+1),mtmp(mpert),xtmp(sq%mx+1,mpert))
             psitmp(:) = sq%xs(0:)
             mtmp = (/(m,m=mlow,mhigh)/)

dcon/makefile

@@ -32,7 +32,8 @@ OBJS = \
 	ode.o \
 	dcon_netcdf.o \
 	resist.o \
-	dcon.o
+	dcon.o \
+	bernstein.o
 
 # targets
 
@@ -84,4 +85,6 @@ free.o: ../equil/global_mod.mod dcon_netcdf.o ode_output.o
 resist.o: dcon_mod.o
 dcon.o: ../equil/equil_mod.mod ../equil/equil_out_mod.mod \
 	../pentrc/inputs.o ../pentrc/pentrc_interface.o \
-	bal.o mercier.o ode.o free.o resist.o
+	bal.o mercier.o ode.o free.o resist.o bernstein.o
+bernstein.o: ../equil/global_mod.mod ../equil/bicube_mod.mod ../equil/spline_mod.mod \
+		dcon_mod.o free.o sing.o 

docs/nametag

@@ -0,0 +1 @@
+snu_offshoot

equil/equil_out.f

@@ -22,6 +22,7 @@
 c-----------------------------------------------------------------------
       MODULE equil_out_mod
       USE global_mod
+
       IMPLICIT NONE
 
       CONTAINS
@@ -878,5 +879,6 @@ SUBROUTINE equil_out_dump
 c     terminate.
 c-----------------------------------------------------------------------
       RETURN
-      END SUBROUTINE equil_out_dump
+      END SUBROUTINE equil_out_dump  
+
       END MODULE equil_out_mod

equil/global.f

@@ -30,6 +30,7 @@ MODULE global_mod
       LOGICAL :: verbose=.TRUE.
       LOGICAL :: use_galgrid=.TRUE.
       LOGICAL :: wv_farwall_flag=.TRUE.
+      LOGICAL :: bernstien_flag=.FALSE.
       REAL(r8) :: psilow=1e-4
       REAL(r8) :: psihigh=1-1e-6
       REAL(r8) :: newq0=0

equil/makefile

@@ -59,4 +59,4 @@ lar.o: inverse.o
 gsec.o: inverse.o
 sol.o: direct.o
 equil.o: read_eq.o lar.o gsec.o sol.o
-equil_out.o: global.o
+equil_out.o: global.o

input/dcon.in

@@ -20,7 +20,7 @@
     mthvac=512           ! Number of points used in splines over poloidal angle at plasma-vacuum interface. Overrides vac.in mth.
     thmax0=1             ! Linear multiplier on the automatic choice of theta integration bounds for high-n ideal ballooning stability computation (strictly, -inf to -inf)
 
-    kin_flag = t         ! Kinetic EL equation (default: false)
+    kin_flag = f         ! Kinetic EL equation (default: false)
     con_flag = t         ! Continue integration through layers (default: false)
     kinfac1 = 1.0        ! Scale factor for energy contribution (default : 1.0)
     kinfac2 = 1.0        ! Scale factor for torque contribution (default : 1.0)
@@ -66,4 +66,5 @@
     bin_bal2=f           ! Binary output for bal_flag functions
 
     netcdf_out=t         ! Replicate ascii dcon.out information in a netcdf file
+    bernstien_flag=t
 /

input/match.in

@@ -5,6 +5,8 @@
 	ideal_flag=t
 	contour_flag=t
 
+	sol_num=1
+
 	res_flag=f
 	ff0=1e6
 	ff1=1e6

match/ideal.f

@@ -305,7 +305,9 @@ END SUBROUTINE ideal_transform
 c-----------------------------------------------------------------------
 c     declarations.
 c-----------------------------------------------------------------------
-      SUBROUTINE ideal_build
+      SUBROUTINE ideal_build(sol_num)
+
+      INTEGER, INTENT(IN) :: sol_num
 
       INTEGER :: istep,ifix,jfix,kfix,ieq,info
       INTEGER, DIMENSION(mpert) :: ipiv
@@ -319,7 +321,8 @@ SUBROUTINE ideal_build
          uedge=0
          uedge(mripple-mlow+1)=1
       ELSE
-         uedge=wt(:,1)
+         WRITE(*,*)"sol_num = ",sol_num
+         uedge=wt(:,sol_num)
       ENDIF
       temp2=soltype(mstep)%u(:,1:mpert,1)
       CALL zgetrf(mpert,mpert,temp2,mpert,ipiv,info)
@@ -357,12 +360,15 @@ END SUBROUTINE ideal_build
 c-----------------------------------------------------------------------
 c     declarations.
 c-----------------------------------------------------------------------
-      SUBROUTINE ideal_write
+      SUBROUTINE ideal_write(sol_num)
+
+      INTEGER, INTENT(IN) :: sol_num
 
       LOGICAL, PARAMETER :: diagnose=.FALSE.
       INTEGER :: ipert,istep,m
       REAL(r8), DIMENSION(mpert) :: singfac
       COMPLEX(r8), DIMENSION(mpert,0:mstep) :: b
+      CHARACTER(128) :: outputname
 c-----------------------------------------------------------------------
 c     compute normal perturbed magnetic field.
 c-----------------------------------------------------------------------
@@ -375,7 +381,14 @@ SUBROUTINE ideal_write
 c     write binary output for graphs.
 c-----------------------------------------------------------------------
       WRITE(*,*)"Write binary output for graphs."
-      CALL bin_open(bin_unit,"solutions.bin","UNKNOWN","REWIND","none")
+
+      IF (sol_num == 1) THEN
+         outputname = "solutions.bin"
+      ELSE
+         WRITE(outputname, '(A,I0,A)') "solutions", sol_num, ".bin"
+      END IF
+
+      CALL bin_open(bin_unit, outputname, "UNKNOWN", "REWIND", "none")
       DO ipert=1,mpert
          DO istep=0,mstep
             WRITE(bin_unit)REAL(psifac(istep),4),REAL(rho(istep),4),

match/match.f

@@ -29,7 +29,7 @@ PROGRAM match_main
      $     big,small,ns,root,scan_flag,root_flag,contour_flag,full,
      $     z_level_type,phi_level,transform_flag,poly_form_type,
      $     poly_test_type,scan0,scan1,nscan,delta0,delta1,condense,
-     $     ripple_flag,mripple
+     $     ripple_flag,mripple,sol_num
 c-----------------------------------------------------------------------
 c     format statements.
 c-----------------------------------------------------------------------
@@ -55,8 +55,8 @@ PROGRAM match_main
 c-----------------------------------------------------------------------
       IF(ideal_flag)THEN
          CALL ideal_transform
-         CALL ideal_build
-         CALL ideal_write
+         CALL ideal_build(sol_num)
+         CALL ideal_write(sol_num)
          CALL ideal_chord
          IF(contour_flag)CALL ideal_contour
       ENDIF

match/match_mod.f

@@ -42,6 +42,7 @@ MODULE match_mod
      $     kin_flag=.FALSE.,con_flag=.FALSE.
       LOGICAL, DIMENSION(:), POINTER :: sing_flag
       INTEGER :: mfix,mhigh,mlow,mpert,mstep,nn,msing,mmatch,mbit,mterm
+      INTEGER :: sol_num = 1
       INTEGER, DIMENSION(:), POINTER :: fixstep
       REAL(r8) :: sfac0=1
       REAL(r8), DIMENSION(:), POINTER :: psifac,rho,q