Jac_calc OMP column version

Makefile.Forthon

@@ -1,13 +1,13 @@
 DEBUG = -v --fargs "-Ofast -DFORTHON -cpp" # -g #--fargs "-C=array" #--fargs "-CB -traceback"
 
 SO = so
-PUBLICHOME = 
+PUBLICHOME =
 BUILDDIR = build
 DIRLIST = com aph api bbb flx grd svr wdf ncl
 MYPYTHON = python3
 FORTHON = Forthon3
 BUILDBASE = -v --build-base $(BUILDDIR)
-INSTALLARGS = --pkgbase uedge $(BUILDBASE) 
+INSTALLARGS = --pkgbase uedge $(BUILDBASE) $(OMPFLAGS)
 
 all: mppl2f90 $(BUILDDIR)/compydep $(BUILDDIR)/grdpydep $(BUILDDIR)/flxpydep $(BUILDDIR)/bbbpydep $(BUILDDIR)/svrpydep $(BUILDDIR)/wdfpydep $(BUILDDIR)/aphpydep $(BUILDDIR)/apipydep $(BUILDDIR)/nclpydep
 	rm -f uedgeC.so
@@ -27,8 +27,8 @@ $(BUILDDIR)/apipydep: $(BUILDDIR)/compydep api/apifcn.F api/apip93.F api/apisorc
 	$(FORTHON) -a $(INSTALLARGS)  -a $(FCOMP) $(DEBUG) --interfacefile api/api.v -d com -f api/apifcn.F api api/apip93.F api/apisorc.F api/fimp.F api/fmombal.F api/inelrates.F api/sputt.F
 	touch $@
 
-$(BUILDDIR)/bbbpydep: $(BUILDDIR)/compydep bbb/boundary.F bbb/convert.F bbb/geometry.F bbb/griddubl.F bbb/oderhs.F bbb/odesetup.F bbb/odesolve.F bbb/potencur.F bbb/turbulence.F bbb/ext_neutrals.F bbb/bbb.v com/com.v 
-	$(FORTHON) -a $(INSTALLARGS)  -a $(FCOMP) --interfacefile bbb/bbb.v $(DEBUG) --macros com/com.v $(READLINE) -d com -f bbb/boundary.F bbb bbb/convert.F bbb/geometry.F bbb/griddubl.F bbb/oderhs.F bbb/odesetup.F bbb/odesolve.F bbb/potencur.F bbb/turbulence.F bbb/ext_neutrals.F 
+$(BUILDDIR)/bbbpydep: $(BUILDDIR)/compydep $(BUILDDIR)/apipydep bbb/boundary.F bbb/convert.F bbb/geometry.F bbb/griddubl.F bbb/oderhs.F bbb/odesetup.F bbb/odesolve.F bbb/potencur.F bbb/turbulence.F bbb/ext_neutrals.F bbb/bbb.v com/com.v
+	$(FORTHON) -a $(INSTALLARGS)  -a $(FCOMP) --interfacefile bbb/bbb.v $(DEBUG) --macros com/com.v $(READLINE) -d com -f bbb/boundary.F bbb bbb/convert.F bbb/geometry.F bbb/griddubl.F bbb/oderhs.F bbb/odesetup.F bbb/odesolve.F bbb/potencur.F bbb/turbulence.F bbb/ext_neutrals.F
 	touch $@
 
 $(BUILDDIR)/flxpydep: $(BUILDDIR)/compydep flx/flxcomp.F flx/flxdriv.F flx/flxread.F flx/flxwrit.F flx/flx.v com/com.v

api/api.v

@@ -157,7 +157,7 @@ work2(nwork2)	_real		# work array for B3VAL
 nwork3	integer			# size of array work3
 work3(nwork3)	_real		# work array for B3INT
 iworki(10)	integer		# work array for B3VAL
-icont		integer	/0/	# input flag for B3VAL
+icont		integer	/0/	+threadprivate # input flag for B3VAL
 kxords_api	integer	/4/	# order of spline fit versus x
 #	kxords_api=4 (default) is cubic interpolation
 kyords_api	integer	/4/	# order of spline fit versus y

api/fimp.m

@@ -95,7 +95,7 @@ c        rate parameters (sigma*v) for ionization, recombination,
       real tmpenonz,tmpinonz,xte,xti,dlogt,fxte,fxti,lrion,lrrec
 c
       real rcxr_zn6, rcxr_zn6b
-      external rcxr_zn6, rcxr_zn6b 
+      external rcxr_zn6, rcxr_zn6b
 c
       rion = 0.
       rrec = 0.
@@ -107,7 +107,7 @@ c        rate parameters (sigma*v) for ionization, recombination,
       xti = log(tmpinonz/ev2)
       dlogt = rtlt(1) - rtlt(0)
 c
-c ... Find index i1 in temperature table such that 
+c ... Find index i1 in temperature table such that
 c                 rtlt(i1) .le. xt .lt. rtlt(i1+1)
 c     or, equivalently,
 c                  rtt(i1) .le. tmp .lt. rtt(i1+1).
@@ -168,13 +168,13 @@ call xerrab("")
 c     Compute rate parameters for transitions from table species ii.
 c
       if (za .lt. zamax) then
-         lrion = 
+         lrion =
      .     (  fy)*((1-fxte)*rtlsa(i1e,j1+1,ii)+fxte*rtlsa(i1e+1,j1+1,ii))
      .   + (1-fy)*((1-fxte)*rtlsa(i1e,j1  ,ii)+fxte*rtlsa(i1e+1,j1  ,ii))
          rion = exp(lrion)
          if (za .eq. 0) return
       endif
-         lrrec = 
+         lrrec =
      .     (  fy)*((1-fxte)*rtlra(i1e,j1+1,ii)+fxte*rtlra(i1e+1,j1+1,ii))
      .   + (1-fy)*((1-fxte)*rtlra(i1e,j1  ,ii)+fxte*rtlra(i1e+1,j1  ,ii))
          rrec = exp(lrrec)
@@ -185,7 +185,7 @@ call xerrab("")
       if ( (iscxfit .gt. 0) .and.
      .     (zn .eq. 6) .and. (za .le. zamax) ) then
          if (iscxfit.ge.1. .and. iscxfit.le.2.) then
-            rcxr = (2.-iscxfit)*rcxr_zn6 (tmpi, za) + 
+            rcxr = (2.-iscxfit)*rcxr_zn6 (tmpi, za) +
      .             (iscxfit-1.)*rcxr_zn6b(tmpi, za)
          endif
 ccc         if (iscxfit .eq. 1) rcxr = rcxr_zn6 (tmpi, za)
@@ -241,8 +241,8 @@ c     Input (neutral hydrogen) temperature, tmp, is in [Joules/AMU].
 c     Output rate parameter (sigma-v) is in [m**3/sec].
 c
 c     This is a modified of the function rcxr_zn6; only za=1 case is
-c     changed to use a (lower) fit guided by plots from A. Pigarov.  
-c     Other za's same as for rxcr_zn6 from thesis by C.F. Maggi (fit 
+c     changed to use a (lower) fit guided by plots from A. Pigarov.
+c     Other za's same as for rxcr_zn6 from thesis by C.F. Maggi (fit
 c     by T. Rognlien)
 c
 c     local variables --
@@ -266,7 +266,7 @@ subroutine readmc(nzdf,mcfilename)
       character*256 mcfilename(*)
       character*256 fname
       Use(Multicharge)
-      Use(Math_problem_size)  # neqmx
+c      Use(Math_problem_size)  # neqmx
       Use(Flags)        # iprint
       Use(Impdata)      #apidir
 
@@ -309,7 +309,7 @@ call xerrab("")
 c     read header --
 *     un*formatted read for header data
       read (nget,'(2a8,i12,4x,a32)') idcod, idtyp, n, id1
-      if (n .lt. 0 .and. iprt_imp_file == 1) then  
+      if (n .lt. 0 .and. iprt_imp_file == 1) then
         if (iprint .ne. 0) then
             write(*,*) '***Impurity file using new 2012 format is ',mcfilename(i)
         endif
@@ -426,7 +426,7 @@ function radmc(zmax, znuc, te, dene, denz, radz)
       real radz(0:zmax)
 c
 c ... Compute the radiation rates, radz(0:zmax), for all charge states
-c     of an impurity with nuclear charge, znuc, and return the total 
+c     of an impurity with nuclear charge, znuc, and return the total
 c     electron energy loss rate, radmc, including both the radiation
 c     and binding energy contributions.
 c
@@ -528,7 +528,7 @@ call xerrab("")
      .  (  fy)*((1-fxt)*rtlqa(i1,j1+1,k0+k)+fxt*rtlqa(i1+1,j1+1,k0+k))
      .+ (1-fy)*((1-fxt)*rtlqa(i1,j1  ,k0+k)+fxt*rtlqa(i1+1,j1  ,k0+k))
          keelz = exp(leelz)
-         fac_rad = 1. 
+         fac_rad = 1.
          if(ispradextrap==1 .and. k==0 .and. te<temintab) then  #extrap below min Te
            fac_rad = (te/(temintab))**6
          endif

api/sputt.m

@@ -5,12 +5,12 @@ SUBROUTINE SYLD96(MATT,MATP,CION,CIZB,CRMB)
 
 c############################################################################
 c ** - modifications by TDR 2/20/98
-c 
+c
 c ** - Coding received from David Elder, Feb. 6, 1998; reduced argument
-c ** - list for subroutine (deleted cneutd,cbombf,cbombz,cebd); deleted 
+c ** - list for subroutine (deleted cneutd,cbombf,cbombz,cebd); deleted
 c ** - params, Use statement for cyield instead of include statement
 c ** - REAL -> real for mppl; ALOG10 -> log10
-c 
+c
 c ** -Input variables:
 c     cion  # integer atomic number of impurity, e.g., for carbon, cion=6
 c     cizb  # integer max charge state of plasma ions; hydrogen cizb=1
@@ -20,7 +20,7 @@ SUBROUTINE SYLD96(MATT,MATP,CION,CIZB,CRMB)
 c     matt  # integer flag giving target material
 c     matp  # integer flag giving plasma material
 c
-c###########################################################################   
+c###########################################################################
 C
 C  *********************************************************************
 C  *                                                                   *
@@ -40,7 +40,7 @@ SUBROUTINE SYLD96(MATT,MATP,CION,CIZB,CRMB)
 cdtr      include    'cyield'
 
       Use(Cyield)   # ceth,cetf,cq,ntars,cidata
-      Use(Math_problem_size) # neqmx
+c      Use(Math_problem_size) # neqmx
       Use(Flags)    # iprint
 
       real ETH(7,12), ETF(7,12), Q(7,12), EBD(12)
@@ -74,7 +74,7 @@ C     THE TWO COMPOUNDS (TITANIUM CARBIDE AND SILICON CARBIDE)
 C     (APPROX 2X).
 C                       LORNE HORTON MAY 93
 C
- 
+
       DATA TARMAT/
      &  ' ALUMINIUM       ',' BERYLLIUM       ',' COPPER          ',
      &  ' GRAPHITE        ',' TITANIUM        ',' IRON            ',
@@ -83,10 +83,10 @@ C     THE TWO COMPOUNDS (TITANIUM CARBIDE AND SILICON CARBIDE)
      &  ' "DEUTERIUM"     ',' "HELIUM"        ',' "NEON"          ',
      &  ' "ARGON"         ',' "OXYGEN"        ',' "CHLORINE"      ',
      &  ' "NITROGEN"      ' /
- 
+
       DATA PLAMAT/
      &  ' H    ',' D    ',' T    ',' HE4  ',' C    ',' SELF ',' O    '/
- 
+
       DATA  ETH/
      &     23.87, 14.86, 12.91, 12.51, 16.32, 24.02, 18.55,
      &     12.2 , 10.0 , 14.69, 13.9 , 28.08, 24.17, 32.71,
@@ -271,17 +271,17 @@ SUBROUTINE SPUTCHEM(IOPTCHEM,E0,TEMP,FLUX,YCHEM)
 c#########################################################################
 c
 c ** -Modified 2/20/98 to use flux in MKS [1/m**2 s] rather than previous
-c ** -[1/cm**2 s]; TDR 
-c 
+c ** -[1/cm**2 s]; TDR
+c
 c ** -Code obtained from David Elder, 2/6/98; originally written by
-c ** -Houyang Guo at JET 
+c ** -Houyang Guo at JET
 c
 c#########################################################################
 C
 C  *********************************************************************
 C  *                                                                   *
 C  *  CHEMICAL SPUTTERING FOR D --> C                                  *
-C  *                                                                   *  
+C  *                                                                   *
 C  *  IOPTCHEM       -  Options for chemical sputtering:               *
 C  *         1       -  Garcia-Rosales' formula (EPS94)                *
 C  *         2       -  according to Pospieszczyk (EPS95)              *
@@ -304,7 +304,7 @@ SUBROUTINE SPUTCHEM(IOPTCHEM,E0,TEMP,FLUX,YCHEM)
 c      Change the flux from MKS to cgs units to use old cgs routines
 c      conversion done on 2/20/98
       FLUX_cgs  = 1e4*FLUX
- 
+
       IF      (IOPTCHEM.EQ.1) THEN
         YCHEM = YGARCIA(E0,TEMP,FLUX_cgs)
       ELSE IF (IOPTCHEM.EQ.2) THEN
@@ -338,7 +338,7 @@ FUNCTION YROTH96(E0,TEMP,FLUX)
 C  *********************************************************************
 C  *                                                                   *
 C  *  CHEMICAL SPUTTERING CALCULATED BY Garcia-Rosales FORMULA        *
-C  *                                                                   *  
+C  *                                                                   *
 C  *  ETHC (eV)  -  Threshold energy for D -> C physical sputtering    *
 C  *  ETFC (eV)  -  Thomas-Fermi energy                                *
 C  *  SNC        -  Stopping power                                     *
@@ -359,14 +359,14 @@ FUNCTION YROTH96(E0,TEMP,FLUX)
 C            Ychem = Ysurf+Ytherm*(1+D*Yphys)
 C ---------------------------------------------------------
 
-C 
+C
 C 1> PHYSICAL SPUTTERING YIELD
 C
       ETHC = 27.0
-      ETFC = 447.0      
+      ETFC = 447.0
       QC   = 0.1
 C
-C  - Stopping Power      
+C  - Stopping Power
 C
       SNC = 0.5*LOG(1.+1.2288*E0/ETFC)/(E0/ETFC
      >    + 0.1728*SQRT(E0/ETFC)
@@ -380,7 +380,7 @@ FUNCTION YROTH96(E0,TEMP,FLUX)
          YPHYS = 0.0
       ENDIF
 C
-C 2> YSURF: Surface Process 
+C 2> YSURF: Surface Process
 C
       CSURF  = 1/(1.+1E13*EXP(-2.45*11604/TEMP))
       CSP3   = CSURF*(2E-32*FLUX+EXP(-1.7*11604/TEMP))
@@ -407,8 +407,8 @@ CW    WRITE(6,*) 'YROTH96 = ',YROTH96
 
       RETURN
       END
-      
-      
+
+
 C -------------
 c ------------------------------------------------------------------------c
       FUNCTION YGARCIA(E0,TEMP,FLUX)
@@ -421,7 +421,7 @@ FUNCTION YGARCIA(E0,TEMP,FLUX)
 C  *********************************************************************
 C  *                                                                   *
 C  *  CHEMICAL SPUTTERING CALCULATED BY Garcia-Rosales FORMULA        *
-C  *                                                                   *  
+C  *                                                                   *
 C  *  ETHC (eV)  -  Threshold energy for D -> C physical sputtering    *
 C  *  ETFC (eV)  -  Thomas-Fermi energy                                *
 C  *  SNC        -  Stopping power                                     *
@@ -432,22 +432,22 @@ FUNCTION YGARCIA(E0,TEMP,FLUX)
 C  *  YCHEM_ATH  -  Athermal mechanism                                 *
 C  *                                                                   *
 C  *********************************************************************
-C            
+C
 
       ETHC = 27.0
-      ETFC = 447.0      
+      ETFC = 447.0
       QC   = 0.1
 C
 C Check for impact energies below threshold
 C
       IF (E0.GT.ETHC) THEN
 C
-C Stopping Power      
+C Stopping Power
 C
       SNC = 0.5*LOG(1.+1.2288*E0/ETFC)/(E0/ETFC
      >    + 0.1728*SQRT(E0/ETFC)
      >    + 0.008*(E0/ETFC)**0.1504)
-C 
+C
 C Physical Sputtering Yield
 C
          YPHYS = QC*SNC*(1-(ETHC/E0)**(2./3.))*(1-ETHC/E0)**2
@@ -471,8 +471,8 @@ CW    WRITE(6,*) 'YTH = ',YCHEM_TH,'YATH = ',YCHEM_ATH
 
       RETURN
       END
-      
-      
+
+
 C -------------
 c ------------------------------------------------------------------------c
 
@@ -484,7 +484,7 @@ FUNCTION YHAASZ(E0,TEMP)
 C  *  - poly. fit: Y = a0 + a1*log(E) + a2*log(E)^2 + a3*log(E)^3      *
 C  *                                                                   *
 C  *********************************************************************
-C            
+C
       IMPLICIT NONE
 
       real    E0,TEMP
@@ -596,7 +596,7 @@ ELSE IF (E0.GT.200.0) THEN
          FITE0 = 200.
       ELSE
          FITE0 = E0
-      ENDIF 
+      ENDIF
 C
       YFIT = 0.0
       DO I = 1,4
@@ -623,7 +623,7 @@ FUNCTION YHAASZ97(E0,TEMP)
 C  *  - poly. fit: Y = a0 + a1*log(E) + a2*log(E)^2 + a3*log(E)^3      *
 C  *                                                                   *
 C  *********************************************************************
-C            
+C
       IMPLICIT NONE
 
       real    E0,TEMP
@@ -724,7 +724,7 @@ ELSE IF (E0.GT.200.0) THEN
          FITE0 = 200.
       ELSE
          FITE0 = E0
-      ENDIF 
+      ENDIF
 C
       YFIT = 0.0
       DO I = 1,4
@@ -755,7 +755,7 @@ FUNCTION YHAASZ97M(E0,TEMP,reducf)
 C  *  Default value for reducf=0.2; change redf_haas                   *
 C  *                                                                   *
 C  *********************************************************************
-C            
+C
       IMPLICIT NONE
 
       real E0, TEMP, reducf
@@ -772,7 +772,7 @@ FUNCTION YHAASZ97M(E0,TEMP,reducf)
          YHAASZ97M = FRAC*YHAASZ97(E0,TEMP)+ (1.-FRAC)*YDAVIS98
       ELSEIF (E0 .LT. 5.) THEN
          YDAVIS98 = reducf/(m2*((TEMP/m1)**2 - 1)**2 + m3)
-         YHAASZ97M = YDAVIS98   
+         YHAASZ97M = YDAVIS98
       ENDIF
 
       RETURN