This file contains the history of edits to the ITG code for 1990 and 1991. The history of more recent edits are in other history files. (This file was original part of the beginning of itg.f, but it got so long that we eventually separated it out.). c This program has evolved from Satoshi Hamaguchi's ETAIMOD program, c which calculates 3D ITG turbulence, as described in the paper c by Hamaguchi and Horton, "Fluctuation Spectrum and Transport from c Ion Temperature Gradient Driven Modes in Sheared Magnetic Fields," c (August, 1989), IFSR-383 and DOE/ET-53088-383. Published in c Physics of Fluids (August, 1990, pp. 1833-1851). c c SEE BELOW SEE BELOW SEE BELOW SEE BELOW SEE BELOW SEE BELOW SEE BELOW c It has been modified to include additional 2cd order FLR corrections c as given in the Eqs. 14-16 of Lee and Tang, Phys. Fluids, (Vol. 31, c 1988, p. 612). These FLR corrections are implemented as Pade c approximations rather than Taylor series approximations, i.e., c as 1/(1+b) rather than as (1-b). c SEE BELOW SEE BELOW SEE BELOW SEE BELOW SEE BELOW SEE BELOW SEE BELOW c c This program evolved from an early version of Wonchull Park's HIB code, c although the convolution section in particular has been rewritten. c c Original version received from Satoshi Hamaguchi on 2-Mar-1990. c c******************************************************************* c c modifications are being made by Bill Dorland (BD), c Greg Hammett (GWH) and Mike Beer (MB): c c MODIFICATION HISTORY: c c BD 1-Oct-1991: Itg.gc is now pseudo-spectral! The subroutine nlps c (NonLinear PseudoSpectral) has replaced the convolutions, and bench- c marked carefully. Still need to install nonlinear dissipative heat c flux and make sure the code vectorizes neatly. Much time is spent c in the FFT calling routines, as the fields are reordered for NAG. c Also, we now solve for T_par directly instead of p_par. c c BD 17-Sep-1991: All radial derivatives now calculated in kx-space. c Need to upgrade to pseudo-spectral (now 1/3 of the way). c c BD 12-Sep-1991: two-temperature update complete. Filters correct. c c BD 23-Aug-1991: q_parallel equation installed nonlinearly. c Filters added to each moment and to the potential of the form c c x -> x/(1 + dt rmu1 b^iexp) c c in order to stabilize high k modes nonlinearly. A more physical c model at high k would be useful. Also, need to redefine thermal c quantities/switches to reflect anisotropic temperature model. c c BD 12-Aug-1991: Changed FLR model. We now approximate several terms c as follows: c c = Gamma_0^(1/2) c = Gamma_0^(1/2) (1 - b/2 (1 - I_1/I_0)) c bar{n} = n - theta(b)*b/2*(1 - I_1/I_0)*T_perp_1 c c where theta(b) = 1 + b/(1 + 3/4*b^2). c c Also have added q_parallel equation linearly. In local limit this c code is very close to results obtained with Mathematica. In slab c limit it compares very well with Linsker's electrostatic cylinder c code, an integral gyrokinetic code detailed in Linsker, Phys. Fluids, c 24, 1485 (1981) obtained from Greg Rewoldt and revived by John c Reynders. c c BD 14-Mar-1991: Consolidated common blocks further, now require all c variables to be explicitly typed to reduce bugs somewhat. Also installed c energy conservation diagnostic, presently guaranteed only for iflr=5. c Benchmarked code against previous runs and got exactly the same results c before I changed the definitions of total energy throughout the code. After c modification, benchmarked again to find no noticeable effect. c c BD 16-Feb-1991: Changed the equations into guiding center coordinates. c c GWH 18-Dec-1990: Added option to do FLR terms exactly using FFT's. c (Although the FFT form was exact in the old physical-space fluid c code based on the Lee-Tang approach, the FFT is only rigorous c to second order in the present second-order gyrofluid code.) c c GWH 14-Dec-1990: Discovered a bug in my implementation of the c Pade approximation that had been used at the APS meeting. The c coefficient of the phi term in the pressure equation had been c (1+etai) too big. Also, I had missed an FLR correction in the c Poisson equation. I'm also trying variations of the Pade c approximations. c c GWH 23-Oct-1990: Added an option to use the additional Lee-Tang c FLR corrections, implemented as Pade approximations rather than c as Taylor series approximations. Hopefully these additional FLR c corrections will improve the comparison with gyrokinetic simulations. c c GWH 23-Oct-1990: I fixed my Landau damping model to have the right c Ti/Te dependence and the proper FLR form. The damping term had c been proportional to (p-phi). It is now proportional to c (p-n) = (p+w*Ti/Te), as it should be. c c BD 19_Sep_1990: There was an error in the error checking routine. c An extra term was kept in the sum of dissipative terms. This is now c removed. Also, the errors are now plotted. This routine must still c be corrected for running with the Hammett-Perkins Landau damping model, c though I have begun this process in this version. c c BD 11_Sep_1990: In order to make movies it is useful to write the data c in the HDF format. The data is now written with the subroutine wframe c in the SDS (scientific data set) format, which allows arrays of data with c greater than two dimensions. c c This routine should be improved to allow more flexibility in choosing the c frames to be written, as the brute force approach seems to generate too c much data: 100 frames of 150 x 200 x 1 data took about 5 megabytes. c c Note: The program is slowed rather dramatically by this option, though I c haven't done a careful timing analysis. Thus, I suggest that movies not c be made as a matter of course (for every run). Until further notice, this c can be accomplished by setting ninterv > nstp. c c Note: A warning is generated during the loading of the hdf library concerning c dynamic common block options used by the C subroutines to generate the c scientific data sets. This should not cause concern. c c BD 10_Sep_1990: The graphics library is now DISSPLA rather than c TV80. The DISSPLA routines are a little slower, but they are portable, c more powerful, and perhaps good enough for immediate public distribution. c The most interesting new feature is the ability to make plots of 2D c surfaces on a 3D set of axes with good perspective. c c BD 26_Aug_1990: The convolution routine was incorrect. It did not c calculate the m=0, n<>0 modes convolutions correctly. The error c was spotted by Wonchull Park and corrected by Park and Dorland. I c checked the new convolution routine by hand on several runs. I c am not sure what the ultimate effect of this change will be, but a c first run shows significant changes in the fluctuation spectrum, c notably in the m=0 mode, and in the value obtained for chi_av. We c should now reexamine the scalings obtained by H&H with extra care. c c Also, the predictor step is now hdt = .5*dt, in order to assure c second-order accuracy as advertised. A benchmark run with moderate c parameters showed little change in the results under this modification. c c GWH: For the model equation df/dt = i omega f, one can show that a c centered two-step method where the first step has hdt=0.5*dt, is c second-order accurate in dt. However, the amplitude of f would slowly c grow due to errors of order dt**4. Using hdt=0.65*dt means the method is c no longer second-order accurate, but it introduces a small amount of c damping of order dt**2. Depending on the size of dt, this damping can c can overcome the growth errors which are of order dt**4. But Bill c saw little change when running with hdt=0.5*dt, so use a centered step. c c BD 23_Aug_1990: Initial conditions given random amplitudes (between c 0. and 1.) and random phases (between 0. and pi). c c BD 22-Aug-1990: Common blocks now referenced via "include" statements. c c Option to choose the parity of the initialization was removed. c (Sine and cosine modes are now chosen.) The magnitudes of these c perturbations should be made random to get away from the very regular and c spiked initial conditions now beng used. c c Option to change time step during run removed. c c Graphics calls concentrated in order to simplify movies and/or a switch to c another graphics library later. c c Finally, several vestigial variables were removed. c c GWH 7-May-1990: Move towards better modularization by replacing c the parameter statements with "include" statements. Common blocks c should be eventually be referenced via "include" statements as well. c c replace "call link" with "open" statements to make more portable. c ITG can now be run on the VAX by splitting up this cosmos file c with csplit.exe, and linking with crayvax.olb, both in c rx15:[itg.greg]. c c GWH 4-Mar-1990: added my fluid model of kinetic effects such c as Landau damping. See Hammett and Perkins, "Fluid Moment Models c for Landau Damping and the Ion Temperature Gradient Instability," c Phys. Rev. Lett. Vol. 64, p. 3019 (June 1990). c c I have only added my Landau-damping model to the equations which c advance the solution in time. There is an energy balance calculation c where it needs to be added as well. c c GWH 3-Mar-1990: minor modifications: added comments, simplified c use of the input namelist file. c c END OF MODIFICATION HISTORY. c*******************************************************************