ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
neowes.F90
Go to the documentation of this file.
1 #ifdef DOIMPURE
2 SUBROUTINE neowes(eq, coreprof, neoclassic, coreimpurity)
3 #else
4 SUBROUTINE neowes(eq, coreprof, neoclassic)
5 #endif
6 
7 !... neoclassical model wrapper
8 !... always use same grid as coreprof
9 !... formulae from Wesson Chapter 4
10 
11  USE itm_constants
12  USE euitm_schemas
13  USE copy_structures
14  USE deallocate_structures
15 
16  IMPLICIT NONE
17 
18  TYPE (type_equilibrium), pointer :: eq(:)
19  TYPE (type_coreprof), pointer :: coreprof(:)
20 #ifdef DOIMPURE
21  TYPE (type_coreimpur), pointer :: coreimpurity(:)
22 #endif
23  TYPE (type_neoclassic), pointer :: neoclassic(:)
24 
25  INTEGER(ITM_I4) :: ion=0,imp=0,icharge=0,iz=0,inucl=0
26  INTEGER(ITM_I4) :: nrho_prof,nion_prof,nrho_eq
27  INTEGER(ITM_I4) :: nrho_imp = 0
28  INTEGER(ITM_I4), SAVE :: &
29  nrho = 0, nion = 0, nimp = 0, ncharge = 0, nspecies = 0
30  REAL(R8) :: a00,b00,r00,rho_tor_max
31  REAL(R8) :: mass,charge
32  INTEGER(ITM_I4), DIMENSION(:), POINTER :: ncharges
33  REAL(R8), DIMENSION(:), POINTER :: &
34  nne,tte,nni,tti,tavg,rlne,rlte,rlni,rlti,qq,zeff,zee, &
35  rhoi,v_i,rnue,rnui,eps12,eps32,mneo
36  REAL(R8), DIMENSION(:), POINTER :: &
37  rho,diffe,diff,chi,vconv,flux,gradn,gradp,eer,sum1,sum2
38  REAL(R8), DIMENSION(:,:), POINTER :: nneo,tneo,zneo
39 
40  REAL(R8) :: kb=itm_ev
41  REAL(R8) :: ee=itm_qe
42  REAL(R8) :: cc=1.0_r8, lcoul=14.0_r8
43  REAL(R8) :: me=itm_me
44  REAL(R8) :: md=itm_md
45 
46  character(len = 132), target :: codename(1) = 'NEOWES'
47  character(len = 132), target :: codeversion(1) = 'CC_Cyprus12'
48 
49 !... find grid size and other stuff from coreprof
50 
51  nrho_prof=SIZE(coreprof(1)%rho_tor)
52  nion_prof=SIZE(coreprof(1)%ni%value, dim=2)
53  nrho_eq=SIZE(eq(1)%profiles_1d%rho_tor)
54 #ifdef DOIMPURE
55  IF (ASSOCIATED(coreimpurity)) THEN
56  IF (SIZE(coreimpurity) > 0) THEN
57  nrho_imp=SIZE(coreimpurity(1)%rho_tor)
58  IF (nrho_imp > 0) THEN
59  nimp=SIZE(coreimpurity(1)%impurity)
60 ! DO imp=1,nimp
61 ! ncharge=SIZE(coreimpurity(1)%impurity(ion)%nz, dim=2)
62 ! END DO
63  END IF
64  END IF
65  END IF
66 #endif
67 
68 !... dimensions and species
69 
70  nrho = nrho_prof
71  nion = nion_prof
72 
73  iz = nion
74 
75 #ifdef DOIMPURE
76  IF (nimp > 0) THEN
77  ALLOCATE(ncharges(nimp))
78 
79  DO imp=1,nimp
80  ncharges(imp) = SIZE(coreimpurity(1)%impurity(imp)%nz, 2)
81  DO icharge=1,ncharges(imp)
82  iz=iz+1
83  END DO
84  END DO
85  END IF
86 #endif
87 
88  nspecies = iz
89 
90 ! WRITE (0,*) 'ions imps species',nion,nimp,iz
91 
92 !... allocations
93 
94  IF (.NOT. ASSOCIATED(neoclassic)) THEN
95  ALLOCATE(neoclassic(1))
96  ALLOCATE(neoclassic(1)%codeparam%codename(1))
97  ALLOCATE(neoclassic(1)%codeparam%codeversion(1))
98  neoclassic(1)%codeparam%codename = codename
99  neoclassic(1)%codeparam%codeversion = codeversion
100 
101  ALLOCATE(neoclassic(1)%rho_tor_norm(nrho))
102  ALLOCATE(neoclassic(1)%rho_tor(nrho))
103  ALLOCATE(neoclassic(1)%sigma(nrho))
104  ALLOCATE(neoclassic(1)%jboot(nrho))
105  ALLOCATE(neoclassic(1)%er(nrho))
106 ! ALLOCATE(neoclassic(1)%fext(nrho,nion,nmoment))
107 ! ALLOCATE(neoclassic(1)%jext(nrho))
108  ALLOCATE(neoclassic(1)%ne_neo%flux(nrho))
109  ALLOCATE(neoclassic(1)%ni_neo%flux(nrho,nion))
110  ALLOCATE(neoclassic(1)%te_neo%flux(nrho))
111  ALLOCATE(neoclassic(1)%ti_neo%flux(nrho,nion))
112  ALLOCATE(neoclassic(1)%ne_neo%diff_eff(nrho))
113  ALLOCATE(neoclassic(1)%ni_neo%diff_eff(nrho,nion))
114  ALLOCATE(neoclassic(1)%te_neo%diff_eff(nrho))
115  ALLOCATE(neoclassic(1)%ti_neo%diff_eff(nrho,nion))
116  ALLOCATE(neoclassic(1)%ne_neo%vconv_eff(nrho))
117  ALLOCATE(neoclassic(1)%ni_neo%vconv_eff(nrho,nion))
118  ALLOCATE(neoclassic(1)%te_neo%vconv_eff(nrho))
119  ALLOCATE(neoclassic(1)%ti_neo%vconv_eff(nrho,nion))
120  ALLOCATE(neoclassic(1)%ti_neo%exchange(nrho,nion))
121  ALLOCATE(neoclassic(1)%ti_neo%qgi(nrho,nion))
122 
123  ALLOCATE(neoclassic(1)%vpol(nrho,nion))
124  ALLOCATE(neoclassic(1)%mtor_neo%flux(nrho))
125  ALLOCATE(neoclassic(1)%mtor_neo%diff_eff(nrho))
126  ALLOCATE(neoclassic(1)%mtor_neo%vconv_eff(nrho))
127 
128  neoclassic(1)%ne_neo%flux=0._r8
129  neoclassic(1)%ni_neo%flux=0._r8
130  neoclassic(1)%te_neo%flux=0._r8
131  neoclassic(1)%ti_neo%flux=0._r8
132  neoclassic(1)%ne_neo%diff_eff=0._r8
133  neoclassic(1)%ni_neo%diff_eff=0._r8
134  neoclassic(1)%te_neo%diff_eff=0._r8
135  neoclassic(1)%ti_neo%diff_eff=0._r8
136  neoclassic(1)%ne_neo%vconv_eff=0._r8
137  neoclassic(1)%ni_neo%vconv_eff=0._r8
138  neoclassic(1)%te_neo%vconv_eff=0._r8
139  neoclassic(1)%ti_neo%vconv_eff=0._r8
140  neoclassic(1)%ti_neo%exchange=0._r8
141  neoclassic(1)%ti_neo%qgi=0._r8
142 
143  neoclassic(1)%vpol=0._r8
144  neoclassic(1)%mtor_neo%flux=0._r8
145  neoclassic(1)%mtor_neo%diff_eff=0._r8
146  neoclassic(1)%mtor_neo%vconv_eff=0._r8
147 
148 #ifdef DOIMPURE
149  IF (nimp > 0) THEN
150  ALLOCATE(neoclassic(1)%nz_neo(nimp))
151  ALLOCATE(neoclassic(1)%tz_neo(nimp))
152  DO imp=1,nimp
153  ncharge=ncharges(imp)
154  ALLOCATE(neoclassic(1)%nz_neo(imp)%flux(nrho,ncharge))
155  ALLOCATE(neoclassic(1)%nz_neo(imp)%diff_eff(nrho,ncharge))
156  ALLOCATE(neoclassic(1)%nz_neo(imp)%vconv_eff(nrho,ncharge))
157  neoclassic(1)%nz_neo(imp)%flux=0._r8
158  neoclassic(1)%nz_neo(imp)%diff_eff=0._r8
159  neoclassic(1)%nz_neo(imp)%vconv_eff=0._r8
160  ALLOCATE(neoclassic(1)%tz_neo(imp)%flux(nrho,ncharge))
161  ALLOCATE(neoclassic(1)%tz_neo(imp)%diff_eff(nrho,ncharge))
162  ALLOCATE(neoclassic(1)%tz_neo(imp)%vconv_eff(nrho,ncharge))
163  ALLOCATE(neoclassic(1)%tz_neo(imp)%exchange(nrho,ncharge))
164  neoclassic(1)%tz_neo(imp)%flux=0._r8
165  neoclassic(1)%tz_neo(imp)%diff_eff=0._r8
166  neoclassic(1)%tz_neo(imp)%vconv_eff=0._r8
167  neoclassic(1)%tz_neo(imp)%exchange=0._r8
168  END DO
169  END IF
170 #endif
171 
172 !... done initialisation
173 
174  END IF
175 
176 !... copy composition over
177 
178  call deallocate_cpo(neoclassic(1)%composition)
179  call copy_cpo(coreprof(1)%composition,neoclassic(1)%composition)
180  call deallocate_cpo(neoclassic(1)%compositions)
181  call copy_cpo(coreprof(1)%compositions,neoclassic(1)%compositions)
182 
183 !... basic parameters
184 
185  a00=eq(1)%eqgeometry%a_minor
186  b00=eq(1)%global_param%toroid_field%b0
187  r00=eq(1)%global_param%toroid_field%r0
188 
189 !... use coreprof grid with equil boundary
190 
191  rho_tor_max=maxval(eq(1)%profiles_1d%rho_tor)
192 
193 !... set neoclassical grid
194 
195  neoclassic(1)%rho_tor=coreprof(1)%rho_tor
196  neoclassic(1)%rho_tor_norm=neoclassic(1)%rho_tor/rho_tor_max
197 
198 !... interpolate profiles
199 
200  rho => neoclassic(1)%rho_tor
201 
202  ALLOCATE(nne(nrho))
203  ALLOCATE(tte(nrho))
204  ALLOCATE(nni(nrho))
205  ALLOCATE(tti(nrho))
206  ALLOCATE(tavg(nrho))
207  ALLOCATE(rlne(nrho))
208  ALLOCATE(rlte(nrho))
209  ALLOCATE(rlni(nrho))
210  ALLOCATE(rlti(nrho))
211  ALLOCATE(qq(nrho))
212  ALLOCATE(zeff(nrho))
213  ALLOCATE(zee(nrho))
214  ALLOCATE(rhoi(nrho))
215  ALLOCATE(v_i(nrho))
216  ALLOCATE(rnue(nrho))
217  ALLOCATE(rnui(nrho))
218  ALLOCATE(eps12(nrho))
219  ALLOCATE(eps32(nrho))
220 
221  ALLOCATE(nneo(nrho,0:nspecies))
222  ALLOCATE(tneo(nrho,0:nspecies))
223  ALLOCATE(zneo(nrho,0:nspecies))
224  ALLOCATE(mneo(0:nspecies))
225 
226  CALL l3interp(coreprof(1)%profiles1d%zeff%value, &
227  coreprof(1)%rho_tor,nrho_prof, zeff,rho,nrho)
228  CALL l3interp(eq(1)%profiles_1d%q,eq(1)%profiles_1d%rho_tor, &
229  nrho_eq, qq,rho,nrho)
230 
231  eps12=sqrt(rho/r00)
232 
233  CALL l3interp(coreprof(1)%ne%value,coreprof(1)%rho_tor, &
234  nrho_prof, nneo(1,0),rho,nrho)
235  CALL l3interp(coreprof(1)%te%value,coreprof(1)%rho_tor, &
236  nrho_prof, tneo(1,0),rho,nrho)
237  zneo(:,0)= - ee
238  mneo(0)=itm_me
239 
240  DO ion=1,nion
241  inucl = coreprof(1)%compositions%ions(ion)%nucindex
242  mass = coreprof(1)%compositions%nuclei(inucl)%amn * itm_amu
243  charge = ee*coreprof(1)%compositions%ions(ion)%zion
244 ! WRITE (0,*) 'ion mass charge',ion,mass,charge
245  CALL l3interp(coreprof(1)%ni%value(1,ion),coreprof(1)%rho_tor, &
246  nrho_prof, nneo(1,ion),rho,nrho)
247  CALL l3interp(coreprof(1)%ti%value(1,ion),coreprof(1)%rho_tor, &
248  nrho_prof, tneo(1,ion),rho,nrho)
249 ! zneo(:,ion)=ee*coreprof(1)%composition%zion(ion)
250 ! mneo(ion)=coreprof(1)%composition%amn(ion) * itm_amu
251  zneo(:,ion)=charge
252  mneo(ion)=mass
253  END DO
254  tavg=sum(nneo(:,1:nion)*tneo(:,1:nion), dim=2)/sum(nneo(:,1:nion), dim=2)
255 
256 #ifdef DOIMPURE
257  IF (nimp > 0) THEN
258  iz = nion
259  DO imp=1,nimp
260  ncharge=ncharges(imp)
261 ! mass=coreimpurity(1)%desc_impur%amn(imp) * itm_amu
262  inucl = coreimpurity(1)%compositions%impurities(imp)%nucindex
263  mass = coreimpurity(1)%compositions%nuclei(inucl)%amn * itm_amu
264 ! WRITE (0,*) 'imp mass',imp,mass
265  DO icharge=1,ncharge
266  iz=iz+1
267  CALL l3interp(coreimpurity(1)%impurity(imp)%nz(1,icharge), &
268  coreimpurity(1)%rho_tor,nrho_imp, &
269  nneo(1,iz),rho,nrho)
270  CALL l3interp(coreimpurity(1)%impurity(imp)%z(1,icharge), &
271  coreimpurity(1)%rho_tor,nrho_imp, &
272  zneo(1,iz),rho,nrho)
273  zneo(:,iz)=ee*zneo(:,iz)
274  tneo(:,iz)=tavg
275  mneo(iz)=mass
276  END DO
277  END DO
278  END IF
279 #endif
280 
281 !... overall species loop
282 
283 ! write (0,*) ' nion nimp nspecies = ',nion,nimp,nspecies
284 ! write (0,*) ' iz imp icharge n T D'
285 
286  DO iz=0,nspecies
287 
288  ion = nion
289 #ifdef DOIMPURE
290  DO imp=1,nimp
291  DO icharge=1,ncharges(imp)
292  ion=ion+1
293  IF (iz == ion) EXIT
294  END DO
295  IF (iz == ion) EXIT
296  END DO
297  IF (iz <= nion) THEN
298  imp=0
299  icharge=0
300  END IF
301 #endif
302 
303  nni=nneo(:,iz)
304  tti=tneo(:,iz)
305 
306 ! WRITE (0,*) 'species min n T',iz,MINVAL(nni),MINVAL(tti)
307 
308  CALL l3deriv(nni,rho,nrho,rlni,rho,nrho)
309  CALL l3deriv(tti,rho,nrho,rlti,rho,nrho)
310  rlni=rlni/nni
311  rlti=rlti/tti
312 
313  mass = mneo(iz)
314  zee = zneo(:,iz)
315 ! rnui=(lcoul/2.09e13_R8)*MAX(nneo(:,0), nni*zee*zee/(ee*ee))/(tti**1.5)
316  IF (imp == 0) THEN
317  rnui=(lcoul/2.09e13_r8)*nni*(zee*zee/(ee*ee))/(tti**1.5)
318  ELSE
319  rnui=(lcoul/2.09e13_r8)* &
320  max( nni*(zee*zee/(ee*ee)), nne )/(tti**1.5)
321  END IF
322 
323  IF (iz == 0) THEN
324  nne=nni
325  tte=tti
326  rlne=rlni
327  rlte=rlti
328  rnue=(lcoul/3.44e11_r8)*zeff*nne/(tte**1.5)
329  rnui=rnue
330  END IF
331 
332 !... derived parameters
333 
334  rhoi=sqrt(cc*cc*mass*kb*tti/(zee*zee*b00*b00))
335  v_i=sqrt(kb*tti/mass)
336 
337 !... electron-ion exchange set to classical
338 !... do conductivity, bootstrap current
339 !... do vpol and electric field at end, save grad p and grad n
340 
341  IF (iz == 0) THEN
342  diff => neoclassic(1)%ne_neo%diff_eff
343  gradn => neoclassic(1)%ne_neo%vconv_eff
344  gradp => neoclassic(1)%ne_neo%flux
345  chi => neoclassic(1)%te_neo%diff_eff
346  flux => neoclassic(1)%te_neo%flux
347 
348  sum1 => neoclassic(1)%er
349  sum2 => neoclassic(1)%te_neo%vconv_eff
350  sum1 = 0.
351  sum2 = 0.
352 
353  neoclassic(1)%sigma = ((1.-eps12)**2.)*nne*ee*ee/(0.51*rnue*itm_me)
354  neoclassic(1)%jboot = nne*kb*tte*(2.44*rlne+0.69*rlte)
355  ELSE IF (iz <= nion) THEN
356  diff => neoclassic(1)%ni_neo%diff_eff(:,iz)
357  gradn => neoclassic(1)%ni_neo%vconv_eff(:,iz)
358  gradp => neoclassic(1)%ni_neo%flux(:,iz)
359  chi => neoclassic(1)%ti_neo%diff_eff(:,iz)
360  flux => neoclassic(1)%ti_neo%flux(:,iz)
361 
362  neoclassic(1)%ti_neo%exchange(:,iz) = &
363  3.*(itm_me/mass)*nne*kb*rnue*(tti-tte)
364  neoclassic(1)%jboot = neoclassic(1)%jboot &
365  + nni*kb*tti*(2.44*rlne-0.42*rlti)
366 #ifdef DOIMPURE
367  ELSE
368  diff => neoclassic(1)%nz_neo(imp)%diff_eff(:,icharge)
369  gradn => neoclassic(1)%nz_neo(imp)%vconv_eff(:,icharge)
370  gradp => neoclassic(1)%nz_neo(imp)%flux(:,icharge)
371  chi => neoclassic(1)%tz_neo(imp)%diff_eff(:,icharge)
372  flux => neoclassic(1)%tz_neo(imp)%flux(:,icharge)
373 #endif
374  END IF
375 
376  gradn = nni*rlni
377  gradp = nni*kb*tti*(rlni+rlti)
378 
379 !... the baseline neoclassical diffusion coefficients
380 !... Pfirsch Schlueter, plateau, banana
381 !... this is greatly simplified, uses the figure from basic lectures
382 !... here we allow addition of a cap on the eps factor
383 !... as near the axis the banana orbits become potato orbits
384 
385  eps32=max(rho,rhoi)/r00
386  eps32=eps32**1.5
387  diff=rhoi*rhoi*rnui*2._r8*qq*qq
388  diff=max(diff, rhoi*rhoi*(v_i/abs(qq*r00))*qq*qq)
389  diff=min(diff, rhoi*rhoi*rnui*qq*qq/eps32)
390 
391 ! write (0,"(5x,3i3,5x,3g12.3)") &
392 ! iz,imp,icharge,nni(nrho/2),tti(nrho/2),diff(nrho/2)
393 
394  chi=3.0*diff
395 
396  flux=nni*kb*tti*(-chi)*rlti
397 
398 !... sums for ambipolarity
399 
400  sum1 = sum1 + zee*diff*gradp/(kb*tti)
401  sum2 = sum2 + nni*zee*zee*diff/(kb*tti)
402 
403  END DO
404 
405 !... finish bootstrap current and electric field
406 
407  eer => neoclassic(1)%er
408 
409  eer = sum1/sum2
410  neoclassic(1)%jboot = - eps12*(qq*r00/((rho+1e-10)*b00))*neoclassic(1)%jboot
411 
412 
413 !... charge correction for diffusion and fluxes in this model
414 !... does it such that charge flux is zero, gives electric field
415 !... D is in diff and gradn is in vconv and grad p is in flux
416 !... correct all species
417 
418  diffe => neoclassic(1)%ne_neo%diff_eff
419 
420  diff => neoclassic(1)%ne_neo%diff_eff
421  vconv => neoclassic(1)%ne_neo%vconv_eff
422  flux => neoclassic(1)%ne_neo%flux
423  gradn => neoclassic(1)%ne_neo%vconv_eff
424  gradp => neoclassic(1)%ne_neo%flux
425 
426  zee = zneo(:,0)
427 
428  flux = -(diff/(kb*tte))*(gradp - nne*zee*eer)
429  vconv = (flux + diff*gradn)/nne
430 
431  DO ion=1,nion
432  diff => neoclassic(1)%ni_neo%diff_eff(:,ion)
433  vconv => neoclassic(1)%ni_neo%vconv_eff(:,ion)
434  flux => neoclassic(1)%ni_neo%flux(:,ion)
435  gradn => neoclassic(1)%ni_neo%vconv_eff(:,ion)
436  gradp => neoclassic(1)%ni_neo%flux(:,ion)
437 
438  zee=zneo(:,ion)
439  nni=nneo(:,ion)
440  tti=tneo(:,ion)
441 
442  neoclassic(1)%vpol(:,ion) = (gradp/(nni*zee) - eer)*cc/b00
443 
444  flux = -(diff/(kb*tti))*(gradp - nni*zee*eer)
445  diff=diffe
446  vconv = (flux + diff*gradn)/nni
447  END DO
448 
449 !... also correct the impurity coefficients
450 
451 #ifdef DOIMPURE
452  iz=nion
453  DO imp=1,nimp
454  DO icharge=1,ncharges(imp)
455  iz=iz+1
456 
457  diff => neoclassic(1)%nz_neo(imp)%diff_eff(:,icharge)
458  vconv => neoclassic(1)%nz_neo(imp)%vconv_eff(:,icharge)
459  flux => neoclassic(1)%nz_neo(imp)%flux(:,icharge)
460  gradn => neoclassic(1)%nz_neo(imp)%vconv_eff(:,icharge)
461  gradp => neoclassic(1)%nz_neo(imp)%flux(:,icharge)
462 
463  zee=zneo(:,iz)
464  nni=nneo(:,iz)
465  tti=tneo(:,iz)
466 
467  flux = -(diff/(kb*tti))*(gradp - nni*zee*eer)
468 ! diff=diffe
469  vconv = (flux + diff*gradn)/nni
470  END DO
471  END DO
472 #endif
473 
474 !... stamp time
475 
476  neoclassic(1)%time=coreprof(1)%time
477 
478 #ifdef DOIMPURE
479  IF (nimp > 0) THEN
480  if(associated(ncharges)) deallocate(ncharges)
481  ENDIF
482 #endif
483  if(associated(nne)) deallocate(nne)
484  if(associated(tte)) deallocate(tte)
485  if(associated(nni)) deallocate(nni)
486  if(associated(tti)) deallocate(tti)
487  if(associated(tavg)) deallocate(tavg)
488  if(associated(rlne)) deallocate(rlne)
489  if(associated(rlte)) deallocate(rlte)
490  if(associated(rlni)) deallocate(rlni)
491  if(associated(rlti)) deallocate(rlti)
492  if(associated(qq)) deallocate(qq)
493  if(associated(zeff)) deallocate(zeff)
494  if(associated(zee)) deallocate(zee)
495  if(associated(rhoi)) deallocate(rhoi)
496  if(associated(v_i)) deallocate(v_i)
497  if(associated(rnue)) deallocate(rnue)
498  if(associated(rnui)) deallocate(rnui)
499  if(associated(eps12)) deallocate(eps12)
500  if(associated(eps32)) deallocate(eps32)
501  if(associated(nneo)) deallocate(nneo)
502  if(associated(tneo)) deallocate(tneo)
503  if(associated(zneo)) deallocate(zneo)
504  if(associated(mneo)) deallocate(mneo)
505 
506 END SUBROUTINE neowes
l3deriv
subroutine l3deriv(y_in, x_in, nr_in, dydx_out, x_out, nr_out)
Definition: l3interp.f90:59
md
subroutine md
Definition: SPAR.f:347
eq
subroutine eq(pcequi, psicon, ncequi, nstep, ngrid,
Definition: Eq_m.f:310
flux
subroutine flux(psitok, rk, zk, nk)
Definition: EQ1_m.f:786
l3interp
subroutine l3interp(y_in, x_in, nr_in, y_out, x_out, nr_out)
Definition: l3interp.f90:1
neowes
subroutine neowes(eq, coreprof, neoclassic)
Definition: neowes.F90:4