diff --git a/Makefile b/Makefile
index 9ffc00138187197da896b0578017acbbdb3d00ce..215e47f28dadaab505b29a145fbb8972028e4649 100644
--- a/Makefile
+++ b/Makefile
@@ -6,14 +6,14 @@ all: dp
# Double precision version
dp: F90 = mpif90
-dp: F90FLAGS = -DNOLAPACK -O3 -xHOST
+dp: F90FLAGS = -O3 -xHOST
dp: EXEC = $(BINDIR)/gyacomo23_dp
dp: dirs src/srcinfo.h
dp: compile
# Single precision version
sp: F90 = mpif90
-sp: F90FLAGS = -DNOLAPACK -DSINGLE_PRECISION -O3 -xHOST
+sp: F90FLAGS = -DSINGLE_PRECISION -O3 -xHOST
sp: EXEC = $(BINDIR)/gyacomo23_sp
sp: dirs src/srcinfo.h
sp: compile
@@ -41,21 +41,21 @@ debug: compile
# For compiling on marconi
marconi: F90 = mpiifort
-marconi: F90FLAGS = -DNOLAPACK -O3 -xHOST
+marconi: F90FLAGS = -O3 -xHOST
marconi: EXEC = $(BINDIR)/gyacomo23_dp
marconi: dirs src/srcinfo.h
marconi: compile
# For compiling on marconi in single prec.
marconi_sp: F90 = mpiifort
-marconi_sp: F90FLAGS = -DNOLAPACK -DSINGLE_PRECISION -O3 -xHOST
+marconi_sp: F90FLAGS = -DSINGLE_PRECISION -O3 -xHOST
marconi_sp: EXEC = $(BINDIR)/gyacomo23_sp
marconi_sp: dirs src/srcinfo.h
marconi_sp: compile
# For compiling on marconi in single prec.
marconi_dbg: F90 = mpiifort
-marconi_dbg: F90FLAGS = -DNOLAPACK -DSINGLE_PRECISION -g -traceback -ftrapuv -warn all -debug all
+marconi_dbg: F90FLAGS = -DSINGLE_PRECISION -g -traceback -ftrapuv -warn all -debug all
marconi_dbg: EXEC = $(BINDIR)/gyacomo23_dbg
marconi_dbg: dirs src/srcinfo.h
marconi_dbg: compile
@@ -93,8 +93,7 @@ gdebug: compile
# test SVD
test_svd: F90 = mpif90
-# test_svd: F90FLAGS = -DTEST_SVD -g -traceback -ftrapuv -warn all -debug all
-test_svd: F90FLAGS = -DTEST_SVD -O3
+test_svd: F90FLAGS = -DLAPACK -DTEST_SVD -O3
test_svd: EXEC = $(BINDIR)/gyacomo23_test_svd
test_svd: dirs src/srcinfo.h
test_svd: compile
@@ -133,7 +132,8 @@ $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/parallel_mod.o
$(OBJDIR)/ppexit.o $(OBJDIR)/prec_const_mod.o \
$(OBJDIR)/processing_mod.o $(OBJDIR)/readinputs.o $(OBJDIR)/restarts_mod.o \
$(OBJDIR)/solve_EM_fields.o $(OBJDIR)/species_mod.o $(OBJDIR)/stepon.o $(OBJDIR)/tesend.o \
-$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
+$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/units_mod.o \
+$(OBJDIR)/CLA_mod.o
# Add finally the definition of the Gyacomo directory
# (helps if the code needs to load a file as for the init Ricci)
@@ -206,7 +206,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/prec_const_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/array_mod.o \
$(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o \
$(OBJDIR)/grid_mod.o $(OBJDIR)/initial_mod.o $(OBJDIR)/model_mod.o \
- $(OBJDIR)/time_integration_mod.o $(OBJDIR)/parallel_mod.o
+ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/parallel_mod.o $(OBJDIR)/units_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnostics_mod.F90 -o $@
$(OBJDIR)/endrun.o : src/endrun.F90 \
@@ -307,7 +307,8 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/readinputs.o : src/readinputs.F90 \
$(OBJDIR)/diagnostics_mod.o $(OBJDIR)/initial_mod.o $(OBJDIR)/model_mod.o \
- $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o \
+ $(OBJDIR)/units_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/readinputs.F90 -o $@
$(OBJDIR)/restarts_mod.o : src/restarts_mod.F90 \
@@ -327,7 +328,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/stepon.o : src/stepon.F90 \
$(OBJDIR)/initial_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field_mod.o \
- $(OBJDIR)/basic_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/grid_mod.o \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/grid_mod.o \
$(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o \
$(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/solve_EM_fields.o\
$(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
@@ -346,6 +347,10 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
+ $(OBJDIR)/units_mod.o : src/units_mod.F90 \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/geometry_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/units_mod.F90 -o $@
+
$(OBJDIR)/CLA_mod.o : src/CLA_mod.F90 \
$(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o \
$(OBJDIR)/time_integration_mod.o
diff --git a/basic_scripts/parameters_CBC.90 b/basic_scripts/parameters_CBC.90
index 288b95aa7d02a6fc1322229f11caf49e9fa6c69b..713e5878889ce49fb834d543eac8699f0497ac0d 100644
--- a/basic_scripts/parameters_CBC.90
+++ b/basic_scripts/parameters_CBC.90
@@ -84,3 +84,13 @@
&TIME_INTEGRATION
numerical_scheme = 'RK4' ! numerical scheme for time integration (RK2,3,4 etc.)
/
+
+&UNITS
+n_ref = 5.0 ! Electron density (x1e19) [m-3]
+T_ref = 2.5 ! Electron temperature [keV]
+R_ref = 1.7 ! Major Radius [m]
+B_ref = 1.5 ! Magnetic field intensity [T]
+m_ref = 1 ! Ion mass [proton mass]
+q_ref = 1 ! Ion charge [elementary charge]
+WRITE_MW = .false. ! Write in std the power output in MW
+/
diff --git a/src/diagnostics_mod.F90 b/src/diagnostics_mod.F90
index cff27d397afc3bf759cb7ca8f57f49703ca7c335..1e8bec5c05c15cfa4fdd6f6a6c061e40bad7f305 100644
--- a/src/diagnostics_mod.F90
+++ b/src/diagnostics_mod.F90
@@ -71,6 +71,7 @@ CONTAINS
USE initial, ONLY: initial_outputinputs
USE time_integration,ONLY: time_integration_outputinputs
USE parallel, ONLY: parallel_outputinputs
+ USE units, ONLY: units_outputinputs
USE futils, ONLY: creatf, creatg, creatd, attach, putfile
IMPLICIT NONE
!input
@@ -114,6 +115,7 @@ CONTAINS
CALL initial_outputinputs(fid)
CALL time_integration_outputinputs(fid)
CALL parallel_outputinputs(fid)
+ CALL units_outputinputs(fid)
! Save STDIN (input file) of this run
IF(jobnum .LE. 99) THEN
WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum
@@ -128,6 +130,7 @@ CONTAINS
USE basic, ONLY: lu_in, chrono_runt, cstep, time, display_h_min_s
USE processing, ONLY: pflux_x, hflux_x
USE parallel, ONLY: my_id
+ USE units, ONLY: WRITE_MW, pow_ref
IMPLICIT NONE
INTEGER, INTENT(in) :: kstep
!! Basic diagnose loop for reading input file, displaying advancement and ending
@@ -150,14 +153,25 @@ CONTAINS
! Terminal info
IF ((kstep .GE. 0) .AND. (MOD(cstep, nsave_0d) == 0) .AND. (my_id .EQ. 0)) THEN
- if(SIZE(pflux_x) .GT. 1) THEN
- WRITE(*,"(A,F8.2,A8,G10.3,A8,G10.3,A8,G10.3,A8,G10.3,A)")&
- '|t = ', time,'| Pxi = ',pflux_x(1),'| Qxi = ',hflux_x(1),'| Pxe = ',pflux_x(2),'| Qxe = ',hflux_x(2),'|'
- else
- WRITE(*,"(A,F8.2,A8,G10.3,A8,G10.3,A)")&
- '|t = ', time,'| Px = ',pflux_x(1),'| Qx = ',hflux_x(1),'|'
- endif
+ IF (WRITE_MW) THEN ! Display MW
+ if(SIZE(pflux_x) .GT. 1) THEN
+ WRITE(*,"(A,F8.2,A8,G10.3,A14,G10.3,A7)")&
+ '|t = ', time,'| Qxi = ',hflux_x(1)*pow_ref,' [MW] | Qxe = ',hflux_x(2)*pow_Ref,' [MW] |'
+ else
+ WRITE(*,"(A,F8.2,A8,G10.3,A8,G10.3,A)")&
+ '|t = ', time,'| Qxi = ',hflux_x(1)*pow_ref,' [MW]|'
+ endif
+ ELSE ! Display normalized code values
+ if(SIZE(pflux_x) .GT. 1) THEN
+ WRITE(*,"(A,F8.2,A8,G10.3,A8,G10.3,A8,G10.3,A8,G10.3,A)")&
+ '|t = ', time,'| Pxi = ',pflux_x(1),'| Qxi = ',hflux_x(1),'| Pxe = ',pflux_x(2),'| Qxe = ',hflux_x(2),'|'
+ else
+ WRITE(*,"(A,F8.2,A8,G10.3,A8,G10.3,A)")&
+ '|t = ', time,'| Pxi = ',pflux_x(1),'| Qxi = ',hflux_x(1),'|'
+ endif
+ ENDIF
ENDIF
+
END SUBROUTINE diagnose
SUBROUTINE diagnose_full(kstep)
diff --git a/src/readinputs.F90 b/src/readinputs.F90
index b6e4bac0ccdd32de5d25238bb336eee68a325734..b5522271179e0f73700780ebcac52ab8a44a7619 100644
--- a/src/readinputs.F90
+++ b/src/readinputs.F90
@@ -1,14 +1,15 @@
SUBROUTINE readinputs
! Additional data specific for a new run
- USE grid, ONLY: grid_readinputs
- USE diagnostics, ONLY: diag_readinputs
+ USE closure, ONLY: closure_readinputs
USE collision, ONLY: collision_readinputs
+ USE diagnostics, ONLY: diag_readinputs
+ USE geometry, ONLY: geometry_readinputs
+ USE grid, ONLY: grid_readinputs
+ USE initial, ONLY: initial_readinputs
USE model, ONLY: model_readinputs
USE species, ONLY: species_readinputs
- USE initial, ONLY: initial_readinputs
USE time_integration, ONLY: time_integration_readinputs
- USE geometry, ONLY: geometry_readinputs
- USE closure, ONLY: closure_readinputs
+ USE units, ONLY: units_readinputs
USE prec_const
IMPLICIT NONE
@@ -44,4 +45,7 @@ SUBROUTINE readinputs
! Load parameters for time integration from input file
CALL time_integration_readinputs
+ ! Load reference profile data in phys. units
+ CALL units_readinputs
+
END SUBROUTINE readinputs
diff --git a/src/units_mod.F90 b/src/units_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..bad4652802ee91d3e30acfd74f22a198d774eb11
--- /dev/null
+++ b/src/units_mod.F90
@@ -0,0 +1,98 @@
+MODULE units
+
+ USE prec_const
+
+ IMPLICIT NONE
+ PRIVATE
+ !! Define the reference values for transforming simulation data into physical units
+ ! The default value, commented here, are corresponding to the #186473 DIII-D discharge at rho=0.95.
+ REAL(xp), PUBLIC, PROTECTED :: n_ref = -1!25 ! density (electrons) [1e19m-3]
+ REAL(xp), PUBLIC, PROTECTED :: T_ref = -1!0.3 ! temperature (electrons) [keV]
+ REAL(xp), PUBLIC, PROTECTED :: R_ref = -1!1.7 ! length (major radius) [m]
+ REAL(xp), PUBLIC, PROTECTED :: B_ref = -1!2.0 ! magnetic field [T]
+ REAL(xp), PUBLIC, PROTECTED :: m_ref = -1!1 ! mass [nucleid mass]
+ REAL(xp), PUBLIC, PROTECTED :: q_ref = -1!1 ! charge [elementary charge]
+ ! Input flag to write heat flux in MW in the std terminal
+ LOGICAL, PUBLIC, PROTECTED :: WRITE_MW = .false.
+ ! Physical constants
+ REAL(xp), PUBLIC, PROTECTED :: kB = 1.4e-23 ! Boltzmann constant [J/K]
+ REAL(xp), PUBLIC, PROTECTED :: mp = 1.7e-27 ! proton mass at rest [kg]
+ REAL(xp), PUBLIC, PROTECTED :: ec = 1.6e-19 ! elementary charge [C]
+
+ ! Derived normalization quantities
+ REAL(xp), PUBLIC, PROTECTED :: c_s0 ! sound speed [m/s]
+ REAL(xp), PUBLIC, PROTECTED :: r_s0 ! sound Larmor radius [m/s]
+ REAL(xp), PUBLIC, PROTECTED :: w_ci ! ion cyclotron frequency [1/s]
+ REAL(xp), PUBLIC, PROTECTED :: A_fs ! flux surface area (torus approx) [m2]
+
+ ! Physical unit converters
+ REAL(xp), PUBLIC, PROTECTED :: hf_ref ! heat flux [W/m2]
+ REAL(xp), PUBLIC, PROTECTED :: pf_ref ! particle flux [1/s/m2]
+ REAL(xp), PUBLIC, PROTECTED :: mf_ref ! momentum flux [kg (m/s)2]
+ REAL(xp), PUBLIC, PROTECTED :: ph_ref ! voltage [V]
+ REAL(xp), PUBLIC, PROTECTED :: pow_ref ! power [MW]
+
+ PUBLIC :: units_readinputs, units_outputinputs
+
+ CONTAINS
+
+ SUBROUTINE units_readinputs
+ ! Read the input parameters
+ USE basic, ONLY: lu_in
+ USE geometry, ONLY: eps, geom
+ USE prec_const, ONLY: xp, pi
+ IMPLICIT NONE
+ REAL(xp) :: T0,n0, p0, m0, q0
+ INTEGER :: fu
+ ! Define namelist to read
+ NAMELIST /UNITS/ n_ref, T_ref, R_ref, B_ref, m_ref, q_ref, WRITE_MW
+ ! Read namelist in input file
+ READ(unit=lu_in,nml=units,iostat=fu)
+ ! Convert in mksa
+ T0 = ec*(1000*T_ref)/kB ! temperature [K]
+ n0 = n_ref*1e19 ! density [m-3]
+ m0 = m_ref*mp ! mass [kg]
+ q0 = q_ref*ec ! charge [C]
+ ! Setup derived normalization quantities
+ c_s0 = sqrt(kB*T0/m0) ! Sound speed [m/s]
+ w_ci = q0*B_ref/m0 ! Ion cyclotron frequ. [Hz]
+ r_s0 = c_s0/w_ci ! Ion Larmor radius [m]
+ SELECT CASE(geom)
+ CASE('z-pinch','Z-pinch','zpinch','Zpinch')
+ A_fs = 2*pi*R_ref*R_ref ! Area of a cylinder of height R_ref [m2]
+ CASE DEFAULT
+ A_fs = 4*pi**2*eps*R_ref**2 ! Area of the toroidal flux surface [m2]
+ END SELECT
+ ! Setup physical unit converters
+ p0 = n0*kB*T0 ! pressure [Pa]
+ hf_ref = p0*c_s0*r_s0**2/R_ref**2 ! Gyrobohm heat flux [W/m2]
+ pow_ref = hf_ref*A_fs/1e6 ! Gyrobohm power [MW]
+ pf_ref = n0*c_s0*r_s0**2/R_ref**2 ! Particle flux [1/s/m2]
+ mf_ref = n0*m0*c_s0**2*r_s0**2/R_ref**2 ! Momentum flux [kg (m/s)2]
+ END SUBROUTINE units_readinputs
+
+
+ SUBROUTINE units_outputinputs(fid)
+ ! Write the input parameters to the results_xx.h5 file
+ USE futils, ONLY: attach, creatd
+ IMPLICIT NONE
+ INTEGER, INTENT(in) :: fid
+ CHARACTER(len=256) :: str
+ WRITE(str,'(a)') '/data/input/units'
+ CALL creatd(fid, 0,(/0/),TRIM(str),'Units Input')
+ CALL attach(fid, TRIM(str), "density (electrons)", n_ref)
+ CALL attach(fid, TRIM(str), "temperature (electrons)", T_ref)
+ CALL attach(fid, TRIM(str), "length (major radius)",R_ref)
+ CALL attach(fid, TRIM(str), "magnetic field (electrons)", B_ref)
+ CALL attach(fid, TRIM(str), "mass (proton)", m_ref)
+ CALL attach(fid, TRIM(str), "charge (elementary)", q_ref)
+ ! Write the mksa converter if units input are provided
+ if ( (hf_ref > 0) .AND. (pf_ref > 0) .AND. (mf_ref > 0)) THEN
+ CALL attach(fid, TRIM(str), "heat flux (gyroBohm)", hf_ref)
+ CALL attach(fid, TRIM(str), "power (gyroBohm)", pow_ref)
+ CALL attach(fid, TRIM(str), "particle flux (gyroBohm)", pf_ref)
+ CALL attach(fid, TRIM(str), "momentum flux (gyroBohm)", mf_ref)
+ ENDIF
+ END SUBROUTINE units_outputinputs
+
+END MODULE units
\ No newline at end of file