text:
HANSARD REVISE * NUMERO 243
Le vendredi 11 juin 1999
LOI SUR LA GESTION DES TERRES DES PREMIERES NATIONS
Motion portant que le debat ne soit plus ajourne
Deuxieme lecture et adoption des amendements du Senat
Deuxieme lecture et adoption des amendements du Senat
LA SEXUALITE ET LA SANTE GENESIQUE
L'ORDRE INTERNATIONAL DES FILLES ET DES FILS DU ROI
LE DECES DE MME SARAH NINGIURUVIK
M. Guy St-Julien
LES CADETS DE L'AIR ET DE L'ARMEE.
LE PARTI REFORMISTE DU CANADA
L'hon. Pierre S. Pettigrew
L'hon. Pierre S. Pettigrew
LA GENDARMERIE ROYALE DU CANADA
M. Gurbax Singh Malhi
L'hon. Pierre S. Pettigrew
L'hon. Arthur C. Eggleton
L'hon. Pierre S. Pettigrew
L'hon. Ronald J. Duhamel
LOI DE 1999 MODIFIANT LES TAXES DE VENTE ET D'ACCISE
Presentation et premiere lecture
LES COMITES DE LA CHAMBRE
HANSARD REVISE * NUMERO 243
Le vendredi 11 juin 1999
La seance est ouverte a 10 heures.
LOI SUR LA GESTION DES TERRES DES PREMIERES NATIONS
MOTION PORTANT QUE LE DEBAT NE SOIT PLUS AJOURNE
L'hon. Don Boudria (leader du gouvernement a la Chambre des communes, Lib.):
Que le debat ne soit plus ajourne.
Plait-il a la Chambre d'adopter la motion?
Que tous ceux qui sont en faveur de la motion veuillent bien dire oui.
Que tous ceux qui sont contre veuillent bien dire non.
A mon avis, les oui l'emportent.
Et plus de cinq deputes s'etant leves:
(La motion, mise aux voix, est adoptee par le vote suivant:)
Je declare la motion adoptee.
Monsieur le President, j'invoque le Reglement.
Je demanderais qu'on retarde la periode des questions de cinq minutes.
A cette periode de l'annee, ma vue baisse un peu.
Est-on d'accord pour proceder ainsi?
DEUXIEME LECTURE ET ADOPTION DES AMENDEMENTS DU SENAT
M. Mike Scott (Skeena, Ref.):
Je trouve cela totalement inacceptable.
La Chambre est-elle prete pour le vote?
Plait-il a la Chambre d'adopter la motion?
(La motion est adoptee et les amendements, lus pour la deuxieme fois, sont adoptes.)
LES TRAVAUX DE LA CHAMBRE
La Chambre a entendu le libelle de la motion.
Lui plait-il d'adopter la motion?
LA SEXUALITE ET LA SANTE GENESIQUE
Mme Beth Phinney (Hamilton Mountain, Lib.):
On estime que, chaque annee, un adolescent sur vingt contracte une maladie transmissible sexuellement.
M. Inky Mark (Dauphin-Swan River, Ref.):
Les Canadiens sont fiers de leur culture.
Le Canada est une societe multiculturelle tres riche et diversifiee.
Nous faisons l'envie du monde entier.
Certains des meilleurs artistes du monde sont originaires du Canada, comme Celine Dion.
L'ORDRE INTERNATIONAL DES FILLES ET DES FILS DU ROI
M. Mac Harb (Ottawa-Centre, Lib.):
Cet organisme a accompli d'innombrables actes charitables au sein de la communaute.
LE DECES DE MME SARAH NINGIURUVIK
M. Guy St-Julien (Abitibi-Baie-James-Nunavik, Lib.):
Elle etait tres connue et impliquee socialement dans sa communaute.
Six compagnies aeriennes n'ont pu se rendre a Kangiqsujuaq en raison de cette reglementation.
M. Bernard Patry (Pierrefonds-Dollard, Lib.):
M. Ken Epp (Elk Island, Ref.):
Et ce qu'ils veulent surtout c'est qu'on reduise les impots.
LES CADETS DE L'AIR ET DE L'ARMEE.
M. Lou Sekora (Port Moody-Coquitlam-Port Coquitlam, Lib.):
Mme Monique Guay (Laurentides, BQ):
J'aimerais feliciter le president Thabo Mbeki, le successeur de Nelson Mandela, pour sa victoire.
Les defis de l'ANC sont encore nombreux et exigeants.
L'hon. Charles Caccia (Davenport, Lib.):
Les mendiants sont plus nombreux.
M. Jay Hill (Prince George-Peace River, Ref.):
Mais hier, un de ces refuges a ete viole.
Les homicides de conjoints sont bien trop frequents dans notre societe.
On compte 100 homicides de conjoints par annee en moyenne au Canada.
J'exhorte tous les deputes a souscrire a cette importante cause.
M. Stan Keyes (Hamilton-Ouest, Lib.):
Monsieur le President, ne vous y trompez pas, nos efforts portent fruit.
Les Canadiens prennent note des avantages qu'apportent nos saines orientations sociales et budgetaires.
Mme Libby Davies (Vancouver-Est, NPD):
Pres de dix ans plus tard, ou en sommes-nous?
La situation des enfants pauvres s'est-elle amelioree?
Il y a cependant un signe d'espoir.
Je compte sur le comite pour qu'il poursuive ses travaux.
Mme Pauline Picard (Drummond, BQ):
Le premier recueil de Roch Napert, Partir est ailleurs , a ete publie en 1997.
Son plus recent recueil, La nuit voyage seule , evoque un univers peuple de riches experiences.
M. Robert Bertrand (Pontiac-Gatineau-Labelle, Lib.):
Nous aurions ete etonnes du contraire.
LE PARTI REFORMISTE DU CANADA
M. Gerald Keddy (South Shore, PC):
Ce genre de comportement est reprehensible et doit cesser a la Chambre.
Voyons le cas du Parti reformiste.
Voila ce qui est reprehensible.
M. Randy White (Langley-Abbotsford, Ref.):
Monsieur le President, avez-vous deja vu une telle peur dans votre vie?
Mme Elinor Caplan (Thornhill, Lib.):
Monsieur le President, dernierement, 13 juifs iraniens ont ete arretes en Iran pour espionnage.
L'arrestation de ces personnes est injustifiee et inacceptable.
M. Randy White (Langley-Abbotsford, Ref.):
Monsieur le President, la periode creuse s'en vient encore une fois a Shawinigan.
L'hon. Herb Gray (vice-premier ministre, Lib.):
Monsieur le President, la question du depute n'est meme pas digne du Parti reformiste.
M. Randy White (Langley-Abbotsford, Ref.):
Le premier ministre n'a pas bien fait les choses dans le passe.
Peut-etre la ministre du Patrimoine pourrait-elle nous aider.
L'hon. Herb Gray (vice-premier ministre, Lib.):
M. Randy White (Langley-Abbotsford, Ref.):
Je demanderais au leader parlementaire de l'opposition de poser directement sa question.
L'hon. Herb Gray (vice-premier ministre, Lib.):
Monsieur le President, premierement, la premisse de la question est inexacte et trompeuse.
Le premier ministre ne possede pas d'actions de la societe en question.
Il les aide, et il le fait de facon tout a fait conforme aux regles.
M. Gurmant Grewal (Surrey-Centre, Ref.):
L'hon. Allan Rock (ministre de la Sante, Lib.):
Ces questions servent a eliminer les risques.
Les consignes de securite appropriees ont ete suivies.
Elle a retrace le sang du donneur.
Elle a fait des tests et les resultats ont ete negatifs.
M. Gurmant Grewal (Surrey-Centre, Ref.):
Malheureusement, quelque chose d'autre a aussi ete prouve.
Le processus de tri du systeme d'approvisionnement en sang est defaillant.
L'hon. Allan Rock (ministre de la Sante, Lib.):
Leur attitude est tout simplement honteuse.
M. Michel Gauthier (Roberval, BQ):
Nous differons d'opinion quant aux consequences de l'entente sur l'union sociale.
M. Michel Gauthier (Roberval, BQ):
Monsieur le President, personne ne peut pretendre serieusement que le Canada est une federation centralisee.
Mais je vais prendre le depute au mot.
Il demande de poser des questions claires aux Quebecois.
M. Stephane Bergeron (Vercheres-Les-Patriotes, BQ):
Quel beau message d'ouverture, monsieur le President.
Le ministre des Affaires intergouvernementales est sensible aux arguments d'autorite.
Il citait pourtant Claude Castonguay, cette semaine, pour defendre l'union sociale.
Ma question s'adresse au ministre des Affaires intergouvernementales.
Non, monsieur le President, ce n'est pas une vision centralisatrice.
Les Quebecois ne sont pas moutonniers.
M. Stephane Bergeron (Vercheres-Les-Patriotes, BQ):
Le ministre des Affaires intergouvernementales a la parole.
Monsieur le President, le Bloc ne fait qu'insulter et lancer des injures.
C'est toi qui lance des injures.
Je veux dire une chose.
Il dit n'importe quoi!
L'hon. Stephane Dion: ...je vais le citer, au sujet du pouvoir federal de depenser:
Je ne suis pas porte a remettre son pouvoir en cause pour deux raisons.
Premierement, le pouvoir de depenser est un pouvoir essentiel des Etats souverains.
M. John Solomon (Regina-Lumsden-Lake Centre, NPD):
Les agriculteurs ne s'interessent par au cirque mediatique.
Ce dont il ont besoin urgent, c'est de l'argent.
M. Joe McGuire (secretaire parlementaire du ministre de l'Agriculture et de l'Agroalimentaire, Lib.):
C'est l'un ou l'autre.
M. John Solomon (Regina-Lumsden-Lake Centre, NPD):
Des inondations sont maintenant a l'origine d'une catastrophe naturelle.
M. Joe McGuire (secretaire parlementaire du ministre de l'Agriculture et de l'Agroalimentaire, Lib.):
M. Peter MacKay (Pictou-Antigonish-Guysborough, PC):
M. Jacques Saada (secretaire parlementaire du solliciteur general du Canada, Lib.):
Ce n'est pas son role.
M. Peter MacKay (Pictou-Antigonish-Guysborough, PC):
L'hon. Herb Gray (vice-premier ministre, Lib.):
Deuxiemement, seuls les frais legaux ont ete payes.
M. Grant McNally (Dewdney-Alouette, Ref.):
Il semble qu'il ait des difficultes d'apprentissage.
Nous avons traite de la meme question, hier soir, a la Chambre.
Ce sont eux qui continuent d'imposer une taxe d'entree aux refugies.
Ils continuent d'imposer une taxe aux refugies, et ils le savent.
En 1996, le Parti liberal a reclame l'abolition de cette taxe d'entree.
Pourquoi le gouvernement est-il determine a arnaquer les refugies?
Je comprends que le Parti reformiste traverse une crise d'identite.
Il ne sait pas ou il s'en va.
M. Bernard Bigras (Rosemont, BQ):
L'honorable ministre du Developpement des ressources humaines a la parole.
L'hon. Pierre S. Pettigrew (ministre du Developpement des ressources humaines, Lib.):
Cela fait un an qu'on le dit.
Ce sont eux, notre priorite, pas ces gens-la.
M. Bernard Bigras (Rosemont, BQ):
L'hon. Pierre S. Pettigrew (ministre du Developpement des ressources humaines, Lib.):
Ils n'ont pas voulu engager un veritable debat sur le fond.
M. Chuck Cadman (Surrey-Nord, Ref.):
M. Jacques Saada (secretaire parlementaire du solliciteur general du Canada, Lib.):
L'honorable depute de Calgary-Est a la parole.
M. Deepak Obhrai (Calgary-Est, Ref.):
Le ministre peut-il me dire ou est ce rapport?
L'hon. Sergio Marchi (ministre du Commerce international, Lib.):
Il sera produit sous peu, monsieur le President.
M. Benoit Sauvageau (Repentigny, BQ):
Ma question s'adresse a l'actuel ministre du Commerce international.
L'hon. Sergio Marchi (ministre du Commerce international, Lib.):
M. Benoit Sauvageau (Repentigny, BQ):
L'hon. Sergio Marchi (ministre du Commerce international, Lib.):
LA GENDARMERIE ROYALE DU CANADA
M. Bill Gilmour (Nanaimo-Alberni, Ref.):
Par contre, la province verse sa juste part.
Quand le solliciteur general va-t-il retablir le financement de base de la GRC?
M. Jacques Saada (secretaire parlementaire du solliciteur general du Canada, Lib.):
Monsieur le President, je pense que mon collegue est un peu en retard.
Cette etude est realisee par une entreprise independante pour nous aider.
M. John Duncan (Ile de Vancouver-Nord, Ref.):
M. Jacques Saada (secretaire parlementaire du solliciteur general du Canada, Lib.):
Ca, le depute n'en parle pas.
L'honorable depute de Sherbrooke a la parole.
M. Serge Cardin (Sherbrooke, BQ):
L'hon. Sergio Marchi (ministre du Commerce international, Lib.):
Dans un dossier d'interet national, les provinces participent, sont consultees et donnent leur avis.
M. Gurbax Singh Malhi (Bramalea-Gore-Malton-Springdale, Lib.):
Mme Brenda Chamberlain (secretaire parlementaire du ministre du Travail, Lib.):
M. Mike Scott (Skeena, Ref.):
L'hon. Jane Stewart (ministre des Affaires indiennes et du Nord canadien, Lib.):
Voila sur quoi nous concentrons nos efforts.
M. Derrek Konrad (Prince Albert, Ref.):
Monsieur le President, ce ministere a le don de laisser les gens en suspens.
Hier, le comite permanent a enfin recu une recommandation.
Un gouvernement liberal anterieur etait present a la conclusion de ce marche boiteux.
L'hon. Jane Stewart (ministre des Affaires indiennes et du Nord canadien, Lib.):
Vous avez manque a votre promesse.
M. Yvon Godin (Acadie-Bathurst, NPD):
L'hon. Pierre S. Pettigrew (ministre du Developpement des ressources humaines, Lib.):
C'est ce qui est au coeur de ce que le premier ministre Bouchard demande.
M. Gordon Earle (Halifax-Ouest, NPD):
Plus de 135 jours plus tard, le ministre n'a toujours pas repondu.
Le ministre est-il reduit au silence par les huiles militaires?
L'hon. Arthur C. Eggleton (ministre de la Defense nationale, Lib.):
Cela a pris un certain temps pour etre traite.
M. Jim Jones (Markham, PC):
Madame la Presidente, ma question s'adresse au vice-premier ministre.
J'ai interjete appel aupres du Commissaire a l'information.
L'hon. Pierre S. Pettigrew (ministre du Developpement des ressources humaines, Lib.):
Monsieur le President, la premisse de la question est absolument fausse.
Le gouvernement respecte la loi du pays concernant l'acces a l'information.
C'est ce que dit la loi et nous la respectons.
M. Jim Jones (Markham, PC):
Le premier ministre se pretend innocent, mais ne fournit aucune preuve corroborante independante.
L'hon. Herb Gray (vice-premier ministre, Lib.):
Cela est tout a fait correct.
Il a agi avec la plus grande integrite.
Il a repondu a toutes les questions clairement, soigneusement et completement.
Mme Sarmite Bulte (Parkdale-High Park, Lib.):
Monsieur le President, ma question s'adresse a la ministre de la Cooperation internationale.
Monsieur le President, il y a beaucoup a faire.
M. Paul Forseth (New Westminster-Coquitlam-Burnaby, Ref.):
C'est pourtant ca le plus important.
Pourquoi expose-t-elle inutilement les enfants aux affres des differends relatifs a leur garde?
Le comite mixte lui a soumis un plan.
Ou est le projet de loi?
M. Michel Bellehumeur (Berthier-Montcalm, BQ):
Ma question est fort simple.
Monsieur le President, je remercie le depute pour sa question.
Cette loi est applaudie par la communaute judiciaire.
Mme Louise Hardy (Yukon, NPD):
Les femmes courent plus de risques d'etre tuees lorsqu'elles quittent leur conjoint.
Veillera-t-elle a ce que les enfants kidnappes soient rendus a leur mere?
Monsieur le President, la violence contre les femmes est un probleme tres grave.
Le gouvernement a pris des mesures pour tenter de le regler.
De toute evidence, le probleme persiste.
La violence faite aux femmes est inacceptable et ne sera pas toleree dans notre societe.
M. Mark Muise (Ouest Nova, PC):
La reponse indiquait que le ministere se preoccupe de la securite des pecheurs.
M. Wayne Easter (secretaire parlementaire du ministre des Peches et des Oceans, Lib.):
Mme Carolyn Bennett (St. Paul's, Lib.):
La recherche est une part decisive pour vaincre le cancer.
Je regrette d'interrompre l'honorable secretaire d'Etat.
J'ai recu avis d'une question de privilege, et je vais l'entendre.
L'hon. Charles Caccia (Davenport, Lib.):
M. John Solomon (Regina-Lumsden-Lake Centre, NPD):
Pourquoi reserve-t-on un traitement particulier a ces deputes, parmi tous les autres?
Mais c'est le cas egalement de bien d'autres deputes.
Pourquoi ce traitement particulier?
L'hon. Don Boudria (leader du gouvernement a la Chambre des communes, Lib.):
Il s'agit du senateur Eugene Whelan.
Je lui souhaite bonne chance dans sa future nomination au Senat.
Une allegation a ete faite et il semble que nous ayons obtenu une explication.
Il ne s'agit pas ici d'un debat.
J'ignore qui l'a fait.
Entre-temps, nous revenons aux affaires curantes.
LOI DE 1999 MODIFIANT LES TAXES DE VENTE ET D'ACCISE
Monsieur le President, j'invoque le Reglement.
Tous les partis sauf un ont donne leur aval.
Je crois savoir que le parti en cause a egalement donne son aval maintenant.
Il n'y a pas de motif inavouable la-dessous.
Je demande le consentement unanime pour adopter la motion suivante:
M. David Pratt (Nepean-Carleton, Lib.):
M. Paul Forseth (New Westminster-Coquitlam-Burnaby, Ref.):
M. Paul Forseth (New Westminster-Coquitlam-Burnaby, Ref.):
M. John Duncan (Ile de Vancouver-Nord, Ref.):
Ces 3 500 signatures ont ete recueillies rapidement, les 6 et 7 juin.
Monsieur le President, je suggere que toutes les questions soient reservees.
Monsieur le President, j'invoque le Reglement.
M. Dick Proctor (Palliser, NPD) propose:
Cette motion en quelque sorte boucle la boucle.
Je voudrais rappeler a la Chambre ce qui s'est passe depuis.
Cete femme s'appelle Lori Foster.
Elle a communique avec le parti en question.
C'etait le Parti reformiste.
On ne l'a pas fait.
Je tiens a recapituler brievement les details.
Je ne veux pas que d'autres personnes soient traitees comme je l'ai ete.
Je partage entierement les sentiments que Mme Foster a exprimes dans sa lettre.
Un parti presente un message s'attaquant a un autre parti.
J'espere que le debat d'aujourd'hui attirera l'attention sur cette question.
C'est un probleme important a Regina.
Ce n'est pas du tout mon intention.
Je pense que nous devons pouvoir citer les faits historiques.
Nous ne pouvons pas recrire l'histoire.
J'attends avec impatience le rapport du comite sur la question.
Le Comite permanent de la procedure et des affaires de la Chambre a declare:
Ce n'est pas la saveur du mois.
Les partis politiques peuvent reagir en quelques heures a des campagnes de publicite les attaquant.
M. Ken Epp (Elk Island, Ref.):
Je partage ce point de vue.
Je crois qu'il est important de lire une fois encore la motion.
Je sais que je devrai rendre des comptes ici.
Je sais que je vais apparaitre de temps a autre a la television.
C'est d'ailleurs pourquoi on ne me voit pas souvent a la television.
Je crois que c'est juste pour tout le monde.
Le depute la-bas a une objection.
Je repete ce que j'ai dit.
C'est l'essence meme d'un debat impartial.
Qu'advient-il si on apporte un amendement?
Qui etablit les limites?
Qu'arrive-t-il s'il s'agit d'un de mes collaborateurs.
Et s'il dit quelque chose?
Et si ca touche ma femme ou ma famille?
La je commence a avoir des reserves.
Je m'oppose contre cette motion, c'est son impraticabilite.
Premierement, elle est trop generale.
Deuxiemement, elle est impraticable.
Qui va fixer les limites?
J'estime plutot que les Canadiens sortent grandis du dialogue.
Ils viennent tres souvent du Quebec, mais aussi de l'est du Canada.
J'exhorte donc de nouveau les deputes a voter contre cette motion.
La derniere chose dont j'aimerais parler porte plus particulierement sur l'affaire Lori Foster.
Elle etait contrariee, car elle croyait manifestement tout ce qu'elle a dit.
Je la respecte pour cela.
Je comprends qu'elle soit demoralisee par l'utilisation de son image.
Toutefois, je ne crois pas qu'elle ait ete mal citee.
C'est une attitude empreinte de jugement et de prejuges.
Je souhaiterais que nous puissions nous en debarrasser.
Nous avons fait cela en toute sincerite et en toute honnetete.
M. Mauril Belanger (secretaire parlementaire de la ministre du Patrimoine canadien, Lib.):
Je donnerai donc mon avis, qui est tout a fait personnel.
Je n'ai jamais ete photogenique.
C'est une raison suffisante pour que j'appuie cette initiative.
Cela aussi devrait peut-etre m'inciter a appuyer cette initiative.
Elle ne dit rien des citations incompletes.
Non seulement c'etait presente hors contexte, mais c'etait incomplet.
Il avait repris la moitie d'une phrase.
La motion du depute devrait traiter de ce genre de choses.
C'est assez restrictif et negatif comme formulation.
Je respecte enormement l'intelligence de l'electorat.
Il n'est pas dupe.
Le depute d'Elk Island a souleve un point interessant.
Il a omis ce point.
Il a omis aussi d'autres points que les collegues vont certainement mentionner.
Representer nos concitoyens a la Chambre des communes est un honneur.
La politique est une profession honorable.
Je pense que nous devrions contre-attaquer.
J'exhorte le depute a suivre une approche positive.
Il aurait alors beaucoup plus de succes et nous serions plus enclins a l'appuyer.
M. Scott Brison (Kings-Hants, PC):
Monsieur le President, je suis heureux d'aborder la motion M-97.
Cela souleve certaines questions en matiere d'ethique.
La Chambre devrait se pencher sur cette question.
A mon avis, la motion est vague.
C'etait la premiere fois que l'on recourait a une publicite aussi franchement negative.
On y voyait le visage de dirigeants politiques quebecois marques d'une croix.
Cette motion devrait etre elargie a d'autres produits aussi.
L'industrie de l'organisation politique a connu une nette croissance.
Je felicite le depute d'avoir presente cette motion.
Je ne vois aucun inconvenient a l'appuyer.
Elle envoie un message utile.
Je crois qu'elle merite notre appui.
M. Stephane Bergeron (Vercheres-Les-Patriotes, BQ):
Deuxiemement, il y a le choix du verbe; on dit s'abstenir.
Cette utilisation avait ete faite par le Parti reformiste.
Cette motion ne fait pas etat de cette realite.
Dans ce cadre, nous pourrions etre d'accord.
Bonnes vacances a tous et toutes, et surtout, bonne periode estivale dans vos circonscriptions.
Monsieur le President, j'invoque le Reglement.
Y a-t-il consentement unanime?
LES COMITES DE LA CHAMBRE
Monsieur le President, j'invoque le Reglement.
Le secretaire parlementaire a-t-il le consentement unanime pour proposer la motion?
L'amendement est-il adopte?
La motion modifiee est-elle adoptee?
(La motion modifiee est adoptee.)
La Chambre reprend l'etude de la motion.
Je souhaite a tous les deputes de belles vacances.
(La seance est levee a 13 h 21.)

ext: f


text:
C    @(#)nrycmp.f	20.4 1/4/99
       subroutine nrycmp
C
C       control Pij on special RZ lines by adjusting Xij.
C $$$
C $$$   Potential improvments could upgrade the sensitivity matrix
C $$$   YCPSEN by the transformation:
C $$$
C $$$   YCPSEN(j+1) = Ti * YCPSEN(j) * Ti,
C $$$
C $$$   where Ti is the elementary diagonal matrix with di = SQRT(ci).
C $$$
C $$$   and ci is the desired ratio dPij/dXij(*) and dPij/dXij(0).
C

      include 'ipfinc/parametr.inc'

      include 'ipfinc/alpha.inc'
      include 'ipfinc/area.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/ecvar.inc'
      include 'ipfinc/intbus.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/prt.inc'
      include 'ipfinc/slnopt.inc'
      include 'ipfinc/trmdbg.inc'
      include 'ipfinc/ycomp.inc'
      include 'ipfinc/ycpsen.inc'
      include 'ipfinc/ycpsln.inc'

        character id * 1
C
C       Use caution in COMPLEX * 16 statement below.  Interface problems
C       may occur if complex arguments are passed as formal parameters
C       to routines which expect COMPLEX * 8, e.g. PIEQIV.
C
        complex * 16 v(4), y(4,4), a(4), s(4), yij, det, yeq(2,2),
     1               oldy11, oldy12, yxx, yxy
        real inew,         imag(MXYCMP),   pij(MXYCMP),  pijadj(MXYCMP),
     1       dx(MXYCMP),   iovld,          vmag(MXYCMP),
     2       dp(MXYCMP),   db(MXYCMP),     serr(MXYCMP), pijvio(MXYCMP),
     3       vmax(MXYCMP), vmin(MXYCMP),   ivio(MXYCMP), iadj(MXYCMP)
        integer status, stat(MXYCMP), statot, debug

        ditot = 0.0
        dytot = 0.0
        kownty = 0
        if (nycomp .eq. 0) go to 900
C            Store line terminal nodal quantities. (They are perturbed
C            during Y-matrix modifications, and distort injection errors
C            for parallel variable compensated lines.)
        statot = 0
        debug = 1
        do 140 jt = 1,nycomp
        stat(jt) = 0
        k1 = kycomp(1,jt)
        k2 = kycomp(2,jt)
        kt = inp2opt(k1)
        mt = inp2opt(k2)
        id = char (kycomp(3,jt))
        ksect = kycomp(4,jt)
        call nrpqv (kt,pk,dpk,qk,dqk,vmag(jt))
        call nrpqv (mt,pk,dpm,qk,dqm,vm)
        serr(jt) = amax1 (abs(dpk),abs(dqk),abs(dpm),abs(dqm))
        dp(jt) = dpk - dpm
C                    Store Vmax, Vmin limits
        vmin(jt) = vlimn(kt)
        vmax(jt) = vlimx(kt)
C                Compute Pij in the following steps:
C                1. Define terminal injections.
        k = 8
        do 100 i = 1,2
        do 100 j = 1,2
           yeq(i,j) = dcmplx (ycomp(k,jt),ycomp(k+1,jt))
           k = k + 2
  100   continue
        v(1) = dcmplx (e(kt),f(kt))
        v(2) = dcmplx (e(mt),f(mt))
        a(1) = yeq(1,1) * v(1) + yeq(1,2) * v(2)
        a(2) = yeq(2,1) * v(1) + yeq(2,2) * v(2)
        xij = ycomp(32,jt)
        bis = ycomp(36,jt)
        dx(jt) = 0.0
        pijadj(jt) = 0.0
C                If this is a passive node sequence, the terminal voltages
C                for section KSECT must be computed by Y-matrix identities.
        if (ksect .eq. 0) then
C                       Store new voltages
           ycomp(44,jt) = dreal (v(1))
           ycomp(45,jt) = dimag (v(1))
           ycomp(46,jt) = dreal (v(2))
           ycomp(47,jt) = dimag (v(2))

           s(1) = v(1) * conjg (a(1))
           pij(jt) = dreal (s(1))

           imag(jt) = dsqrt((dreal(a(1)))**2 + (dimag(a(1)))**2)

           imag(jt) = sign (imag(jt),pij(jt))

        else
C
C       1. Initialize the 4-node Y-matrix. The rows and columns are
C          ordered 1 (Node 1), 2 (Node 2), 3 (Node x), and 4 (Node y).
C
           do 110 j = 1,4
           do 110 i = 1,4
              y(i,j) = 0.0
  110      continue
C
C       2. Append subsystem (1-x), (x-y), and (y-2).
C
           k = 16
           do 120 i = 1,3,2
           do 120 j = 1,3,2
              y(i,j) = y(i,j) + dcmplx (ycomp(k,jt),ycomp(k+1,jt))
              k = k + 2
  120      continue

           yij = dcmplx (-1.0,0.0)/dcmplx(0.0,xij)
           y(3,4) = y(3,4) + yij
           y(4,3) = y(4,3) + yij
           y(3,3) = y(3,3) + dcmplx (0.0,bis) - yij
           y(4,4) = -yij

           k = 24
           do 130 i = 2,4,2
           do 130 j = 2,4,2
           y(i,j) = y(i,j) + dcmplx (ycomp(k,jt),ycomp(k+1,jt))
  130      k = k + 2
C
C       3. Modify row 3 if nodes 1-x are coalesced.
C
           if (y(1,1) .eq. 0.0) then
              a(3) = a(1)
           else
              a(3) = -y(3,1) * v(1)
           endif
C
C       4. Modify row 4 if nodes 2-y are coalesced.
C
           if (y(2,2) .eq. 0.0) then
              a(4) = a(2)
           else
              a(4) = -y(4,2) * v(2)
           endif
C
C          Compute voltages
C
           det = y(3,3) * y(4,4) - y(3,4) * y(4,3)
           v(3) = (y(4,4) * a(3) - y(3,4) * a(4)) / det
           v(4) = (-y(4,3) * a(3) + y(3,3) * a(4)) / det
C
C          Store new voltages
C
           ycomp(44,jt) = dreal (v(3))
           ycomp(45,jt) = dimag (v(3))
           ycomp(46,jt) = dreal (v(4))
           ycomp(47,jt) = dimag (v(4))

           yxy = dcmplx (-1.0,0.0)/dcmplx(0.0,xij)
           yxx = dcmplx (0.0,bis) - yxy

           a(3) = yxx * v(3) + yxy * v(4)
           s(3) = v(3) * conjg (a(3))
           pij(jt) = dreal (s(3))

           imag(jt) = dsqrt((dreal(a(3)))**2 + (dimag(a(3)))**2)
           imag(jt) = sign (imag(jt),pij(jt))

        endif
C
C       Compute residuals
C
        dpij = dim (pij(jt),ycomp(6,jt)) - dim (ycomp(5,jt),pij(jt))
        dytot = dytot + abs (dpij)

        if (ycomp(40,jt) .gt. 0.0) then
           ditot = ditot + dim (imag(jt),ycomp(40,jt)) +
     1                     dim (-ycomp(40,jt),imag(jt))
        endif
C
C       If injection error dominates, ignore X and dPij/dXij updating
C       this iteration.
C
        if (abs (dpij) .gt. option(6)) then
           da = dp(jt)/dpij
           if ( abs (da) .gt. 2.0) stat(jt) = stat(jt) + 100
        endif

        if (abs (serr(jt)) .gt. option(5)) stat(jt) = stat(jt) + 200
C
C       YCPSCR(1,*) = Pij(o)
C       YCPSCR(2,*) = Pin(j-1)
C       YCPSCR(3,*) = dXij(j-1)
C       YCPSCR(4,*) = dPij(0)/dXij
C       YCPSCR(5,*) = dPij(*)/dXij
C       YCPSCR(6,*) = Pij (best) - Pmax
C
        dpdxol = ycomp(42,jt)
        dpdxne = ycomp(42,jt)

        if (ycpscr(4,jt) .eq. 0.0) then
           ycpscr(1,jt) = pij(jt)
           ycpscr(2,jt) = pij(jt)
           ycpscr(3,jt) = 0.0
           ycpscr(4,jt) = ycomp(42,jt)
           ycpscr(5,jt) = ycomp(42,jt)
           ycpscr(6,jt) = pij(jt) - ycomp(5,jt)
        else if (ycpscr(3,jt) .eq. 0.0) then
           ycpscr(2,jt) = pij(jt)
           ycpscr(6,jt) = pij(jt) - ycomp(5,jt)
        else if (ycpscr(3,jt) .ne. 0.0) then
           dpdxne = (pij(jt) - ycpscr(2,jt)) / ycpscr(3,jt)
C
C          Conditionally update dPij/dXij(*) if two conditions are
C          met:
C
C          1. Pij(j) is closer to scheduled Pmax.
C          2. The ratio of the measured dPij/dXij is "close" to the
C             anticipated value dPij/dXij(*).
C
           ratio = dpdxne / ycpscr(5,jt)
           if (ratio .lt. 0.25 .or .ratio .gt. 4.0) then
              ycpscr(2,jt) = pij(jt)
              ycpscr(3,jt) = 0.0
              ycpscr(6,jt) = pij(jt) - ycomp(5,jt)
           else
              if (abs (pij(jt) - ycomp(5,jt)) .le. abs (ycpscr(6,jt)))
     1           then
                 ycpscr(5,jt) = ycpscr(4,jt) +
     1                         (ycomp(5,jt) - ycpscr(1,jt)) *
     2                         (dpdxne - ycpscr(4,jt)) /
     3                         (pij(jt) - ycpscr(1,jt))
                 ycpscr(6,jt) = pij(jt) - ycomp(5,jt)
              endif
           endif
        endif
C
C       Exterpolate sensitivity dPij/dXij
C
        dpdxop = ycpscr(4,jt) +
     1           (pij(jt) - ycpscr(1,jt)) *
     2           (ycpscr(4,jt) - ycpscr(5,jt)) /
     3           (ycpscr(1,jt) - ycomp(5,jt))
C
C       Set flag to abandon DX, DB adjustment if DPDXOP is significantly
C       different from DPDXNE.
C
        ratio = dpdxne / dpdxop
        if (ratio .lt. 0.25 .or .ratio .gt. 4.0) then
           stat(jt) = stat(jt) + 500
        endif
        if (ycpscr(3,jt) .ne. 0 .and. stat(jt) .lt. 100) then
           ycomp(42,jt) = dpdxop
        endif

  140   statot = statot + stat(jt)
C
C       Conditionally modify sensitivities to conform with DPDXOP.
C
        if (statot .eq. 0 .and. statot .ne. 0) then
C
C       Logic disabled 4 Oct 1985 - WLP
C
           if (debug .ne. 0) then
              do 146 jt = 1,nycomp
              if (jt .eq. 1) then
                 write (dbug,142) jt,(j,ycpsen(jt,j),j=1,nycomp)
  142            format (' NRYCMP/YCPSEN (old) - row ',i3,
     1                   5(1x,i4,1x,e12.4) /
     2              (31x,5(1x,i4,1x,e12.4)))
                 if (iterm .ne. 0) then
                    write (*,142) jt,(j,ycpsen(jt,j),j=1,nycomp)
                 endif
              else
                 write (dbug,144) jt,(j,ycpsen(jt,j),j=1,nycomp)
  144            format ('                       row ',i3,
     1                   5(1x,i4,1x,e12.4) /
     2              (31x,5(1x,i4,1x,e12.4)))
                 if (iterm .ne. 0) then
                    write (*,144) jt,(j,ycpsen(jt,j),j=1,nycomp)
                 endif
              endif
  146         continue
           endif
           do 160 jt = 1, nycomp
           ratio = sqrt (ycomp(42,jt) / ycpsen(jt,jt))
           do 150 j = 1, nycomp
           if (j .ne. jt) then
              ycpsen(jt,j) = ycpsen(jt,j) * ratio
              ycpsen(j,jt) = ycpsen(j,jt) * ratio
           endif
  150      continue
  160      continue
           do 170 jt = 1, nycomp
           ycpsen(jt,jt) = ycomp(42,jt)
  170      continue
           if (debug .ne. 0) then
              do 176 jt = 1,nycomp
              if (jt .eq. 1) then
                 write (dbug,172) jt,(j,ycpsen(jt,j),j=1,nycomp)
  172            format (' NRYCMP/YCPSEN (new) - row ',i3,
     1                   5(1x,i4,1x,e12.4) /
     2              (31x,5(1x,i4,1x,e12.4)))
                 if (iterm .ne. 0) then
                    write (*,172) jt,(j,ycpsen(jt,j),j=1,nycomp)
                 endif
              else
                 write (dbug,174) jt,(j,ycpsen(jt,j),j=1,nycomp)
  174            format ('                       row ',i3,
     1                   5(1x,i4,1x,e12.4) /
     2              (31x,5(1x,i4,1x,e12.4)))
                 if (iterm .ne. 0) then
                    write (*,174) jt,(j,ycpsen(jt,j),j=1,nycomp)
                 endif
              endif
  176         continue
           endif
        endif
C       Obtain DX,DB using LP: Minimize Z,
C
C       Z = CP * (P_VIOLATION) + CV * (V_VIOLATION) + CI * (I_VIOLATION)
C
        call ycplp (pij,pijadj,imag,vmag,vmin,vmax,dx,db,pijvio,
     1              iadj,ivio,z,status)

        do 300 jt = 1,nycomp

        k1 = kycomp(1,jt)
        k2 = kycomp(2,jt)

        kt = inp2opt(k1)
        mt = inp2opt(k2)

        xij = ycomp(32,jt)
        bis = ycomp(36,jt)

        xijnew = xij + dx(jt)
        bisnew = bis + db(jt)
        pijnew = pij(jt) + pijadj(jt)

        k = 8
        do 190 i = 1,2
        do 190 j = 1,2
        yeq(i,j) = dcmplx (ycomp(k,jt),ycomp(k+1,jt))
  190   k = k + 2
C
C       Compute residuals
C
        inew = imag(jt) + iadj(jt)
        if (ycomp(40,jt) .gt. 0.0) then
           iovld = dim (abs(inew),ycomp(40,jt))
        else
           iovld = 0.0
        endif

        if (idswb .ne. 0) then
           write (dbug,200) jt, status, stat(jt), kt, mt, z,
     1        pij(jt), pijadj(jt), pijnew, pijvio(jt),
     2        xij, dx(jt), xijnew, ycomp(34,jt), ycomp(35,jt),
     3        bis, db(jt), bisnew, ycomp(38,jt), ycomp(39,jt),
     4        imag(jt), iadj(jt), inew, ivio(jt), iovld,
     5        dpdxol, dpdxne, dpdxop
  200      format (' NRYCMP JT,STATUS,STAT,KT,MT,Z       ',5i5,5x,e10.3
     1           / '        PIJ,DPIJ,PIJNEW,PIJVIO       ',4f10.5
     2           / '        XIJ,DX,XIJNEW,XIJMIN,XIJMAX  ',5f10.5
     3           / '        BIS,DB,BISNEW,BISMIN,BISMAX, ',5f10.5
     4           / '        IOLD,DI,INEW,IVIO,IOVLD      ',5f10.5
     5           / '        dP/dX(0),(j-1),(j)           ',3e10.3 )
           if (iterm .ne. 0) then
              write (*,200) jt, status, stat(jt), kt, mt, z,
     1           pij(jt), pijadj(jt), pijnew, pijvio(jt),
     2           xij, dx(jt), xijnew, ycomp(34,jt), ycomp(35,jt),
     3           bis, db(jt), bisnew, ycomp(38,jt), ycomp(39,jt),
     4           imag(jt), iadj(jt), inew, ivio(jt), iovld,
     5           dpdxol, dpdxne, dpdxop

           endif
        endif

        if (abs (dx(jt)) .le. 1.0e-5 .or. statot .ge. 100) then

        else

           kownty = kownty + 1
C
C          Store new sensitivity values:
C
C          YCPSCR(2,*) = Pij(j-1)
C          YCPSCR(3,*) = dXij(j-1)
C
           ycpscr(2,jt) = pij(jt)
           ycpscr(3,jt) = dx(jt)
C
C          Store new admittances
C
           ycomp(32,jt) = xijnew
           ycomp(36,jt) = bisnew
C
C          Modify Y-matrix
C
C       1. Initialize the 4-node Y-matrix. The rows and columns are
C          ordered 1 (Node 1), 2 (Node 2), 3 (Node x), and 4 (Node y).
C
           do 210 i = 1,4
           do 210 j = 1,4
  210      y(i,j) = 0.0
C
C       2. Append subsystem (1-x), (x-y), and (y-2).
C
           k = 16
           do 220 i = 1,3,2
           do 220 j = 1,3,2
           y(i,j) = y(i,j) + dcmplx (ycomp(k,jt),ycomp(k+1,jt))
  220      k = k + 2

           if (xijnew .eq. 0.0) xijnew = 1.0e-6
           yij = dcmplx (-1.0,0.0)/dcmplx(0.0,xijnew)
           y(3,4) = y(3,4) + yij
           y(4,3) = y(4,3) + yij
           y(3,3) = y(3,3) + dcmplx (0.0,bisnew) - yij
           y(4,4) = -yij

           k = 24
           do 230 i = 2,4,2
           do 230 j = 2,4,2
           y(i,j) = y(i,j) + dcmplx (ycomp(k,jt),ycomp(k+1,jt))
  230      k = k + 2
C
C       3. If Y(1,1) = 0, then coalesce nodes x and 1;
C          otherwise eliminate node x.
C
           if (y(1,1) .eq. 0.0) then
C
C             r1' = r1 + r3;
C             c1' = c1 + c3
C
C             Y(1,1) = Y(1,1) + Y(1,3) + Y(3,1) + Y(3,3)
C             Y(1,4) = Y(1,4) + Y(3,4)
C             Y(4,1) = Y(4,1) + Y(4,3)
C
              y(1,1) = y(3,3)
              y(1,4) = y(3,4)
              y(4,1) = y(4,3)
           else
              y(3,1) = y(3,1) / y(3,3)
              y(3,4) = y(3,4) / y(3,3)
              y(4,4) = y(4,4) - y(4,3) * y(3,4)
              y(4,1) = -y(4,3) * y(3,1)
              y(1,1) = y(1,1) - y(1,3) * y(3,1)
              y(1,4) = -y(1,3) * y(3,4)
           endif
C
C       4. If Y(2,2) = 0, then coalesce nodes y and 2;
C          otherwise eliminate node x.
C
           if (y(2,2) .eq. 0.0) then
C
C             r2' = r2 + r4;
C             c2' = c2 + c4
C
C             Y(2,2) = Y(2,2) + Y(2,4) + Y(4,2) + Y(4,4)
C             Y(1,2) = Y(1,2) + Y(1,4)
C             Y(2,1) = Y(2,1) + Y(4,1)
C
              y(2,2) = y(4,4)
              y(1,2) = y(1,4)
              y(2,1) = y(4,1)
           else
              y(4,1) = y(4,1) / y(4,4)
              y(4,2) = y(4,2) / y(4,4)
              y(1,1) = y(1,1) - y(1,4) * y(4,1)
              y(1,2) = y(1,2) - y(1,4) * y(4,2)
              y(2,1) = y(2,1) - y(2,4) * y(4,1)
              y(2,2) = y(2,2) - y(2,4) * y(4,2)
           endif
C
C          Y(1,1):Y(2,2) now contains the new Y-matrix values.
C
           k = 8
           do 232 i = 1,2
           do 232 j = 1,2
           ycomp(k,jt) = dreal (y(i,j))
           ycomp(k+1,jt) = dimag (y(i,j))
  232      k = k + 2
           do 240 l = km(kt), km(kt)-1+kmlen(kt)
           if (ikmu(l) .eq. mt) then
              oldy12 = dcmplx (gkmu(l),bkmu(l))
              l1 = l
              gkmu(l) = gkmu(l) - dreal (yeq(1,2)) + dreal (y(1,2))
              bkmu(l) = bkmu(l) - dimag (yeq(1,2)) + dimag (y(1,2))
       endif
  240  continue
              oldy11 = dcmplx (gkku(kt) ,bkku(kt))
              l2 = l
              gkku(kt) = gkku(kt) - dreal (yeq(1,1)) + dreal (y(1,1))
              bkku(kt) = bkku(kt) - dimag (yeq(1,1)) + dimag (y(1,1))
           if (idswb .ne. 0) then
               write (dbug,250) kt,mt,oldy11,oldy12,gkku(kt),bkmu(kt),
     1           gkmu(l1),bkmu(l1)
  250         format (' Y-matrix adjustment : Old Y11 ',2i6,2e12.5,
     1              /,'                       Old Y12 ',12x,2e12.5,
     2              /,'                       New Y11 ',12x,2e12.5,
     3              /,'                       New Y12 ',12x,2e12.5)
           endif

           do 260 l = km(mt), km(mt)-1+kmlen(mt)
           if (ikmu(l) .eq. kt) then
              oldy12 = dcmplx (gkmu(l), bkmu(l))
              l1 = l
              gkmu(l) = gkmu(l) - dreal (yeq(2,1)) + dreal (y(2,1))
              bkmu(l) = bkmu(l) - dimag (yeq(2,1)) + dimag (y(2,1))
           endif
  260      continue
              oldy11 = dcmplx (gkku(mt),bkku(mt))      ! diagonal
              gkku(mt) = gkku(mt) - dreal (yeq(2,2)) + dreal (y(2,2))
              bkku(mt) = bkku(mt) - dimag (yeq(2,2)) + dimag (y(2,2))
           if (idswb .ne. 0) then
              write (dbug,250) mt,kt,oldy11,oldy12,gkku(mt),bkku(mt),
     1           gkmu(l1), bkku(l1)
           endif
C                    Modify TIE admittance if K1 - K2 is also a tie line
           if (kycomp(7,jt) .ne. 0) then
              iax = kycomp(7,jt)
              if (tie(1,iax) .eq. kt .and. tie(7,iax) .eq. mt) then

                 oldy11 = dcmplx (tie(3,iax),tie(4,iax))
                 oldy12 = dcmplx (tie(5,iax),tie(6,iax))

                 tie(5,iax) = tie(5,iax) - dreal (yeq(1,2))
     1                                   + dreal (y(1,2))
                 tie(6,iax) = tie(6,iax) - dimag (yeq(1,2))
     2                                   + dimag (y(1,2))
                 tie(3,iax) = tie(3,iax) - dreal (yeq(1,1))
     1                                   + dreal (y(1,1))
                 tie(4,iax) = tie(4,iax) - dimag (yeq(1,1))
     1                                   + dimag (y(1,1))
              else if (tie(1,iax) .eq. mt .and. tie(7,iax) .eq. kt)
     1           then

                 oldy11 = dcmplx (tie(3,iax),tie(4,iax))
                 oldy12 = dcmplx (tie(5,iax),tie(6,iax))

                 tie(5,iax) = tie(5,iax) - dreal (yeq(2,1))
     1                                   + dreal (y(2,1))
                 tie(6,iax) = tie(6,iax) - dimag (yeq(2,1))
     2                                   + dimag (y(2,1))
                 tie(3,iax) = tie(3,iax) - dreal (yeq(2,2))
     1                                   + dreal (y(2,2))
                 tie(4,iax) = tie(4,iax) - dimag (yeq(2,2))
     1                                   + dimag (y(2,2))
              else
                 call erexit
              endif
              if (idswb .ne. 0) then
                 write (dbug,270) iax,ifix(sngl(tie(1,iax))),
     &              ifix(sngl(tie(7,iax))),
     1              oldy11,oldy12,tie(3,iax),tie(4,iax),
     2              tie(5,iax),tie(6,iax)
  270            format (' Tie adjustment',i6,': Old Y11 ',2i6,2e12.5,
     1              /,'                          Old Y12 ',12x,2e12.5,
     2              /,'                          New Y11 ',12x,2e12.5,
     3              /,'                          New Y12 ',12x,2e12.5)
              endif
           endif
        endif
  300   continue
        write (outbuf,310) ditot, dytot, kownty
  310   format (t53,2f13.5,i7)
        call prtout(1)
  900   continue
        return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
C     LAST UPDATE 21/02/94
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE STOREC
      IMPLICIT NONE
C
C Purpose: Converts coordinates from standard to reciprocal space
C
C Calls   0:
C Called by: CTILT
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C FIX include file:
C
      INCLUDE 'FIXPAR.COM'
C
C Local variables:
C
      INTEGER I
      REAL CMU,CHI,TEMP,R,Z,D
C
C-----------------------------------------------------------------------
      DO 10 I=1,NPTS
         CMU = COS(ATAN2(SXY(1,I),SDD))
         CHI = ATAN2(SXY(2,I),SQRT(SXY(1,I)**2+SDD**2))
         TEMP = 1.0 - CMU*COS(CHI)
         IF(TEMP.LT.0.0)TEMP = 0.0
         D = SQRT(2.0*TEMP)/WAVE
         Z = (SIN(TILT)*TEMP+COS(TILT)*SIN(CHI))/WAVE
         R = D*D - Z*Z
         IF(R.LT.0.0)R = 0.0
         R = SQRT(R)
         R = SIGN(R,SXY(1,I))
         RZ(1,I) = R
         RZ(2,I) = Z
 10   CONTINUE
      RETURN
C
      END


module lowflow

  implicit none

  private

  public :: get_lowflow_terms, finish_lowflow_terms
  public :: phineo, dphidth, mach_lab

  real, dimension (:), allocatable :: dphidth
  real, dimension (:,:,:,:,:), allocatable :: coefs
  real, dimension (:,:), allocatable :: phineo, dcoefsdr
  real, dimension (:,:,:), allocatable :: dcoefsdth
  real, dimension (:), allocatable :: rad_neo, theta_neo
  real, dimension (:,:), allocatable :: dens_neo, temp_neo
  real, dimension (:), allocatable :: mach_lab

contains
  subroutine finish_lowflow_terms
    implicit none
    !<DD>May actually be able to deallocate most of these much
    !    earlier as most only used during initialisation
    if(allocated(dphidth)) deallocate(dphidth)
    if(allocated(coefs)) deallocate(coefs)
    if(allocated(phineo)) deallocate(phineo)
    if(allocated(dcoefsdr)) deallocate(dcoefsdr, dcoefsdth)
    if(allocated(rad_neo)) deallocate(rad_neo)
    if(allocated(theta_neo)) deallocate(theta_neo)
    if(allocated(mach_lab)) deallocate(mach_lab)
    if(allocated(dens_neo)) deallocate(dens_neo)
    if(allocated(temp_neo)) deallocate(temp_neo)
  end subroutine finish_lowflow_terms

  subroutine get_lowflow_terms (theta, al, energy, bmag, dHdEc, dHdxic, vpadHdEc, dHdrc, &
       dHdthc, hneoc, dphidrc, dphidthc, phi_neo, lf_default, lf_decompose)

    use mp, only: proc0, broadcast
    use le_grids, only: w, wl
    use theta_grid, only: delthet, jacob, ntgrid, Rplot
    use file_utils, only: open_output_file, close_output_file, get_unused_unit
    use geometry, only: rmaj
    use species, only: spec

    implicit none

    real, dimension (:), intent (in) :: theta, bmag
    real, dimension (:), intent (in) ::  al
    real, dimension (:), intent (in) :: energy
    real, dimension (:,:,:,:,:), intent (out) :: dHdec, dHdxic, dHdrc, dHdthc, vpadHdEc, hneoc
    real, dimension (:), intent (out) :: dphidthc, dphidrc, phi_neo
    logical, intent (in) :: lf_default, lf_decompose

    real :: rgeo2, drlt, drln

    real, dimension (:,:,:,:,:), allocatable :: hneo
    real, dimension (:,:,:,:), allocatable :: dHdxi, dHdE, vpadHdE, dHdr, dHdth
    real, dimension (:,:,:), allocatable :: legp
    real, dimension (:,:), allocatable :: xi, emax, dphidr
    real, dimension (:,:,:,:), allocatable :: chebyp1, chebyp2
    real, dimension (:), allocatable :: dl_over_b, drad_neo
    real, dimension (:,:), allocatable :: transport
    real, dimension (:), allocatable :: phineo2
    real, dimension (:,:), allocatable :: upar2, qpar2
    real, dimension (:,:,:), allocatable :: upar, qpar, afac, bfac
    real, dimension (:,:), allocatable :: pflx, qflx, vflx, qparB, upar1, uparB, upar_over_B
    real, dimension (:), allocatable ::  radius, jboot, phi2, vtor, upar0
    real, dimension (size(energy)) :: w_r

    integer :: il, ie, is, ns, nc, nl, nr, ig, ixi, ir, ir_loc
    integer :: ntheta, nlambda, nenergy, nxi
    integer, save :: neo_unit, neot_unit

!    logical, save :: initialized = .false.
    logical, dimension (:,:), allocatable :: forbid

!    if (initialized) return
!    initialized = .true.

    ntheta = size(theta)
    nlambda = size(al)
    nenergy = size(energy)
    nxi = 2*nlambda-1

    if (.not. allocated(xi)) allocate (xi(ntheta,nxi))
    if (.not. allocated(forbid)) allocate(forbid(ntheta,nxi))
    forbid = .false. ; xi = 0.0

    ! get xi=vpar/v from lambda=(mu/E)*Bnorm
    do il = 1, nlambda-1
       do ig = 1, ntheta
          if (1.0-al(il)*bmag(ig) < 0.0) then
             forbid(ig,il) = .true.
             forbid(ig,2*nlambda-il) = .true.
          else
             xi(ig,il) = -sqrt(1.0-al(il)*bmag(ig))
             xi(ig,2*nlambda-il) = -xi(ig,il)
          end if
       end do
    end do

    ! ns is nspec, nc is nenergy, nl is nlambda,
    ! nr is number of radii
    ! taken from neo
    call read_neocoefs (theta, ns, nc, nl, nr, ir_loc)

    if (.not. allocated(chebyp1)) allocate (chebyp1(nr,nenergy,0:nc-1,ns), chebyp2(nr,nenergy,0:nc-1,ns))
    if (.not. allocated(legp)) allocate (legp(ntheta,nxi,0:nl+1))
    if (.not. allocated(hneo)) allocate (hneo(nr,ntheta,nxi,nenergy,ns))
    if (.not. allocated(dHdr)) allocate (   dHdr(ntheta,nxi,nenergy,ns))
    if (.not. allocated(dHdth)) allocate (  dHdth(ntheta,nxi,nenergy,ns))
    if (.not. allocated(dHdxi)) allocate (  dHdxi(ntheta,nxi,nenergy,ns))
    if (.not. allocated(dHdE)) allocate (   dHdE(ntheta,nxi,nenergy,ns))
    if (.not. allocated(vpadHdE)) allocate (vpadHdE(ntheta,nxi,nenergy,ns))
    if (.not. allocated (dphidr)) allocate (dphidr(ntheta,ns))
    if (.not. allocated (dl_over_b)) allocate (dl_over_b(ntheta))
    if (.not. allocated (emax)) allocate (emax(nr,ns))

    ! better to be taken from neo
    if (proc0) write (*,*) '# make sure ENERGY_MAX=16.0 in NEO INPUT file'
    ! emax is ENERGY_MAX (from NEO input file) times v_{ts}^2

    ! two ways to deal with radial variation of temperature in
    ! NEO energy variable: 1st (lf_default) is to take it into
    ! account when constructing distribution function, i.e. H(EGS2,r);
    ! 2nd is to construct H(ENEO) and deal with it later by
    ! including a temp gradient term in wstarfac in dist_fn.fpp
    if (lf_default) then
       ! v_{ts} is a function of radius, so we need to convert emax
       ! to its equivalent value using v_{ts} from the center radius.
       ! this is necessary because we will be taking radial derivatives of
       ! the distribution function with v fixed, not v/v_t(r) fixed.
       do is = 1, ns
          emax(2,is) = 16.0 ! this is EMAX for center grid point
          emax(1,is) = emax(2,is)*temp_neo(1,is)
          emax(3,is) = emax(2,is)*temp_neo(3,is)
!          emax(1,is) = emax(2,is)*(1.0-spec(is)%tprim*(rad_neo(1)-rad_neo(2)))
!          emax(3,is) = emax(2,is)*(1.0-spec(is)%tprim*(rad_neo(3)-rad_neo(2)))
       end do
    else
       emax = 16.0
    end if

    allocate (drad_neo(nr)) ; drad_neo = 0.
    drad_neo = rad_neo - rad_neo(2)

    ! get legendre polynomials on gs2 pitch-angle grid
    legp = 0.0
    do ixi = 1, nxi
       do ig = 1, ntheta
          if (.not. forbid(ig,ixi)) call legendre (xi(ig,ixi), legp(ig,ixi,:))
       end do
    end do

    ! get chebyshev polynomials of the first and second kinds on the gs2 energy grid
    ! note the argument of the chebyshev polynomials is z = 2 * sqrt(energy/emax) - 1
    do is = 1, ns
       do ie = 1, nenergy
          do ir = 1, nr
             call chebyshev (zfnc(energy(ie),emax(ir,is)), chebyp1(ir,ie,:,is), 1)
             call chebyshev (zfnc(energy(ie),emax(ir,is)), chebyp2(ir,ie,:,is), 2)
          end do
       end do
    end do

    ! first, get hneo = F_1/F_0 = H_1/F_0 - Z_s * e * Phi / T_s
    ! H_1/F_0 is constructed from Legendre and Chebyshev polynomials
    do is = 1, ns
       do ig = 1, ntheta
          do ir = 1, nr
             call get_H (coefs(ir,ig,:,:,is), legp(ig,:,:), chebyp1(ir,:,:,is), &
                  hneo(ir,ig,:,:,is), phineo(ir,ig), ig, is)
          end do
       end do
    end do

    ! calculate and write to file some useful quantities like Upar and Qpar
    dl_over_b = delthet*jacob
    dl_over_b = dl_over_b/sum(dl_over_b)

    if (.not. allocated(pflx)) allocate (pflx(nr,ns), qflx(nr,ns), vflx(nr,ns), upar1(nr,ns))
    if (.not. allocated(upar)) then
       allocate (phineo2(nr))
       allocate (upar2(nr,ns), qpar2(nr,ns))
       allocate (uparB(nr,ns), qparB(nr,ns))
       allocate (upar(nr,ns,ntheta), qpar(nr,ns,ntheta))
       allocate (upar_over_b(nr,ns))
       upar = 0. ; qpar = 0.
       upar2 = 0. ; qpar2 = 0. ; phineo2 = 0.
       uparB = 0. ; qparB = 0.
       upar_over_B = 0.
    end if
    if (.not. allocated(transport)) allocate (transport(nr,(5+ns*8)))  !JPL, NEO output for all radius
    if (.not. allocated(radius)) allocate (radius(nr), phi2(nr), jboot(nr), vtor(nr), upar0(nr))
    if (.not. allocated(mach_lab)) allocate (mach_lab(ns))

    if (proc0) then
       call get_unused_unit (neo_unit)
       ! read the diagnosis in NEO from "neo_transport.out"
       open (unit=neo_unit, file='neo_transport.out', status='old', action='read')
       do ir=1, nr
          read (neo_unit,*) transport(ir,:)
          radius(ir) = transport(ir,1)
          phi2(ir) = transport(ir,2)
          jboot(ir) = transport(ir,3)
          vtor(ir) = transport(ir,4)
          upar0(ir) = transport(ir,5)
          do is = 1, ns
             pflx(ir,is) = transport(ir,5+(is-1)*8+1)/sqrt(2.0)  ! NEO uses vt=sqrt(T/m) while GS2 uses vt=sqrt(2T/m)
             qflx(ir,is) = transport(ir,5+(is-1)*8+2)/(sqrt(2.0)**3)
             vflx(ir,is) = transport(ir,5+(is-1)*8+3)/(sqrt(2.0)**2)
             upar1(ir,is) = transport(ir,5+(is-1)*8+4)/sqrt(2.0)
          end do
       end do
       close (neo_unit)

       ! evaluate the neoclassical parallel heat flux and parallel velcotiy in GS2 geometry
       do ir=1, nr
          do is=1, ns
             do ie = 1, nenergy
                w_r(ie) = w(ie) * exp(energy(ie)*(1.0-1.0/temp_neo(ir,is))) / temp_neo(ir,is)**1.5
                do ixi = 1, nxi
                   il = min(ixi, nxi+1-ixi)
                   do ig = 1, ntheta
                      upar(ir,is,ig) = upar(ir,is,ig) + hneo(ir,ig,ixi,ie,is) &
                           * sqrt(energy(ie))*xi(ig,ixi) * w_r(ie) * wl(-ntgrid+ig-1,il)
                      qpar(ir,is,ig) = qpar(ir,is,ig) + hneo(ir,ig,ixi,ie,is) &
                           * sqrt(energy(ie))*xi(ig,ixi)*(energy(ie)-2.5*temp_neo(ir,is)) &
                           * dens_neo(ir,is) * w_r(ie) * wl(-ntgrid+ig-1,il)
                   end do
                end do
             end do
!              drlt = drad_neo(ir)*spec(is)%tprim ; drln = drad_neo(ir)*spec(is)%fprim
!              do ie = 1, nenergy
!                 w_r(ie) = w(ie) * exp(energy(ie)*(1.0-1.0/(1.0-drlt))) / (1.0-drlt)**1.5
!                 do ixi = 1, nxi
!                    il = min(ixi, nxi+1-ixi)
!                    do ig = 1, ntheta
!                       upar(ir,is,ig) = upar(ir,is,ig) + hneo(ir,ig,ixi,ie,is) &
!                            * sqrt(energy(ie))*xi(ig,ixi) * w_r(ie) * wl(-ntgrid+ig-1,il)
!                       qpar(ir,is,ig) = qpar(ir,is,ig) + hneo(ir,ig,ixi,ie,is) &
!                            * sqrt(energy(ie))*xi(ig,ixi)*(energy(ie)-2.5*(1.0-drlt)) &
!                            * (1.0-drln) * w_r(ie) * wl(-ntgrid+ig-1,il)
!                    end do
!                 end do
!              end do
          end do
       end do

       do is = 1, ns
          do ir = 1, nr
             uparB(ir,is) = sum(upar(ir,is,:)*dl_over_b*bmag)
             qparB(ir,is) = sum(qpar(ir,is,:)*dl_over_b*bmag)
             upar_over_B(ir,is) = sum(upar(ir,is,:)*dl_over_b/bmag)
             upar2(ir,is) = sum(upar(ir,is,:)**2*dl_over_b)
             qpar2(ir,is) = sum(qpar(ir,is,:)**2*dl_over_b)
             phineo2(ir) = sum(phineo(ir,:)**2*dl_over_b)
          end do
       end do

       rgeo2 = 0.
       do ig = 1, ntheta
          rgeo2 = rgeo2 + dl_over_b(ig)*Rplot(-ntgrid+ig-1)**2
       end do

       ! note that for total toroidal angular momentum (including diamagnetic and ExB) to be zero,
       ! upar_over_b/rgeo2 = -(omega_phi * a / v_{t,ref}) * (a/R0)

       call open_output_file (neot_unit,".neotransp")
      ! print out in ".neotransp" file
       write (neot_unit,fmt='(a14,a8,12a14)') "# 1) rad", "2) spec", "3) pflx", "4) qflx", "5) vflx", &
            "6) qparflx", "7) uparB(GS2)", "8) upar1(NEO)", "9) <phi**2>", "10) bootstrap", "11) upar_o_B", &
            "12) upar2", "13) qpar2", "14) phi2"
       do ir=1, nr
          do is = 1, ns
             write (neot_unit,fmt='(e14.5,i8,12e14.5)') radius(ir), is, pflx(ir,is), qflx(ir,is), vflx(ir,is), &
                  qparB(ir,is),uparB(ir,is), upar1(ir,is)*sqrt(temp_neo(ir,is)), &
                  phi2(ir), jboot(ir), rmaj*upar_over_B(ir,is)/rgeo2, upar2(ir,is), qpar2(ir,is), phineo2(ir)
             if (ir==2)  mach_lab(is)=rmaj*upar_over_B(ir,is)/rgeo2
          end do
       end do
       call close_output_file (neot_unit)
    end if

    call broadcast(mach_lab)

    do ir = 1, nr
       do is = 1, ns
          call broadcast (upar(ir,is,:))
          call broadcast (qpar(ir,is,:))
       end do
    end do

    allocate (afac(nr,ns,ntheta), bfac(nr,ns,ntheta))

    ! Here only to set up F_1/F_0 = xi*(v/vth)*(a + b*(v/vth)^2)
    if (lf_decompose) then

       ! for testing
!       upar = 0.
!       qpar(1,:,:) = qpar(2,:,:) ; qpar(3,:,:) = qpar(2,:,:)
!       qpar = 0.
       ! added after doing the qpar=0 test (to decouple upar and dupar/dr effects)
!       upar(1,:,:) = upar(2,:,:) ; upar(3,:,:) = upar(2,:,:)
       ! added after doing the qpar=0,du/dr=0 test...am eliminating the du/dtheta effect
!       do is = 1, ns
!          do ir = 1, nr
!             upar(ir,is,:) = sum(upar(ir,is,:)*dl_over_b)
!          end do
!       end do

       do is = 1, ns
          do ir = 1, nr
             drlt = drad_neo(ir)*spec(is)%tprim ; drln = drad_neo(ir)*spec(is)%fprim
             afac(ir,is,:) = 2.*(upar(ir,is,:)-qpar(ir,is,:)/(1.0-drln) - drlt*upar(ir,is,:))/(1.0-drlt)**2
             bfac(ir,is,:) = 0.8*qpar(ir,is,:) / ((1.0-drln) * (1.0-drlt)**3)
!             afac(ir,is,:) = 2.*(upar(ir,is,:)-qpar(ir,is,:)/(1.0-drln)/(1.0-drlt))/sqrt(1.0-drlt)
!             bfac(ir,is,:) = 0.8*qpar(ir,is,:) / ((1.0-drln) * (1.0-drlt)**1.5)
          end do
       end do

       coefs = 0.0
       do ig = 1, ntheta
          coefs(:,ig,1,0,:) = sqrt(emax)*(0.5*afac(:,:,ig) + 0.3125*bfac(:,:,ig)*emax)
          coefs(:,ig,1,1,:) = 2.*sqrt(emax)*(0.25*afac(:,:,ig) + 0.234375*bfac(:,:,ig)*emax)
          coefs(:,ig,1,2,:) = 2.*0.09375*emax**1.5*bfac(:,:,ig)
          coefs(:,ig,1,3,:) = 2.*0.015625*emax**1.5*bfac(:,:,ig)
       end do

       ! first, get hneo = F_1/F_0 = H_1/F_0 - Z_s * e * Phi / T_s
       ! H_1/F_0 is constructed from Legendre and Chebyshev polynomials
       do is = 1, ns
          do ig = 1, ntheta
             do ir = 1, nr
                call get_H (coefs(ir,ig,:,:,is), legp(ig,:,:), chebyp1(ir,:,:,is), &
                     hneo(ir,ig,:,:,is), phineo(ir,ig), ig, is)
             end do
          end do
       end do

       upar = 0. ; qpar = 0.
       do ir = 1, nr
          do is = 1, ns
             drlt = drad_neo(ir)*spec(is)%tprim ; drln = drad_neo(ir)*spec(is)%fprim
             do ie = 1, nenergy
                w_r(ie) = w(ie) * exp(energy(ie)*(1.0-1.0/(1.0-drlt))) / (1.0-drlt)**1.5
                do ixi = 1, nxi
                   il = min(ixi, nxi+1-ixi)
                   do ig = 1, ntheta
                      upar(ir,is,ig) = upar(ir,is,ig) + hneo(ir,ig,ixi,ie,is) &
                           * sqrt(energy(ie))*xi(ig,ixi) * w_r(ie) * wl(-ntgrid+ig-1,il)
                      qpar(ir,is,ig) = qpar(ir,is,ig) + hneo(ir,ig,ixi,ie,is) &
                           * sqrt(energy(ie))*xi(ig,ixi)*(energy(ie)-2.5*(1.0-drlt)) &
                           * (1.0-drln) * w_r(ie) * wl(-ntgrid+ig-1,il)
                   end do
                end do
             end do
          end do
       end do

       do is = 1, ns
          do ir = 1, nr
             upar2(ir,is) = sum(upar(ir,is,:)**2*dl_over_b)
             qpar2(ir,is) = sum(qpar(ir,is,:)**2*dl_over_b)
          end do
       end do

    end if

    ! get dphi/dr at GS2 radius (center of 3 NEO radii)
    ! note that radgrad takes derivative of e*phi/Tref,
    ! not the derivative of Z_s*e*phi/T_s
    do ig = 1, ntheta
       ! Note that we are neglecting the term prop. to dTref/dr, which is fine
       ! as long as Tref in NEO is the same for all radii.
       ! This will require a changing TEMP with radius if the species temperature
       ! changes with radius.
       call get_radgrad (phineo(:,ig), rad_neo, ir_loc, dphidr(ig,1))
    end do

    ! get radial derivative of F_1/F_0 at fixed (xi,E)
    ! note that F_1/F_0 = H_1/F_0 - Z_s * e * Phi / T_s
!     do is = 1, ns
!        do ig = 1, ntheta
!           do ie = 0, nc-1
!              do ixi = 0, nl
!                 ! get radial derivative of spectral coefficients of H_1/H_0
!                 call get_radgrad (coefs(:,ig,ixi,ie,is), rad_neo, ir_loc, dcoefsdr(ixi,ie))
!              end do
!           end do
!           call get_gradH (dcoefsdr, dphidr(ig,is), legp(ig,:,:), chebyp1(ir_loc,:,:,is), dHdr(ig,:,:,is))
!        end do
!     end do

    ! get theta derivative of F_1/F_0 at fixed (xi,E)
    do is = 1, ns
       do ie = 0, nc-1
          do ixi = 0, nl
             ! get theta derivative of spectral coefficients of F_1/F_0
             call get_thgrad (coefs(ir_loc,:,ixi,ie,is), theta, dcoefsdth(:,ixi,ie))
          end do
       end do
       ! note that dphidth is calculated in read_neocoefs
       do ig = 1, ntheta
          call get_gradH (dcoefsdth(ig,:,:), dphidth(ig)*spec(is)%zt, legp(ig,:,:), chebyp1(ir_loc,:,:,is), dHdth(ig,:,:,is))
       end do
    end do

    do is = 1, ns
       do ig = 1, ntheta
          ! get dH/dxi at fixed E and dH/dE at fixed xi
          call get_dHdxi (coefs(ir_loc,ig,:,:,is), legp(ig,:,:), chebyp1(ir_loc,:,:,is), xi(ig,:), dHdxi(ig,:,:,is))
          call get_dHdE (coefs(ir_loc,ig,:,:,is), legp(ig,:,:), chebyp1(ir_loc,:,:,is), chebyp2(ir_loc,:,:,is), &
               energy(:), emax(ir_loc,is), dHdE(ig,:,:,is))
       end do
    end do

    ! get radial derivative of F_1/F_0 at fixed (xi,E)
    ! note that F_1/F_0 = H_1/F_0 - Z_s * e * Phi / T_s
    do is = 1, ns
       do ig = 1, ntheta
!          do ie = 0, nc-1
!             do ixi = 0, nl
!                ! get radial derivative of spectral coefficients of H_1/H_0
!                call get_radgrad (coefs(:,ig,ixi,ie,is), rad_neo, ir_loc, dcoefsdr(ixi,ie))
          do ie = 1, nenergy
             do ixi = 1, nxi
                call get_radgrad (hneo(:,ig,ixi,ie,is), rad_neo, ir_loc, dHdr(ig,ixi,ie,is))
             end do
          end do
!          call get_gradH (dcoefsdr, dphidr(ig,is), legp(ig,:,:), chebyp1(ir_loc,:,:,is), dHdr(ig,:,:,is))
       end do
    end do

    ! get dH/dtheta and dH/dr
!     do is = 1, ns
!        do ie = 1, nenergy
!           do ixi = 1, nxi
!              do ig = 1, ntheta
!                 call get_dHdr (hneo(:,ig,ixi,ie,is), rad_neo, ir_loc, dHdr(ig,ixi,ie,is))
!              end do
!              call get_dHdth (hneo(ir_loc,:,ixi,ie,is), theta, dHdth(:,ixi,ie,is))
!           end do
!        end do
!     end do

! deprecated ->    ! this is actually dphi/dr - phi*d(ln T)/dr, where phi = e*Phi/Tref
!    dphidrc(1:ntheta-1) = 0.5*(dphidr(1:ntheta-1,1) + dphidr(2:ntheta,1))/spec(1)%zt
    dphidrc(1:ntheta-1) = 0.5*(dphidr(1:ntheta-1,1) + dphidr(2:ntheta,1))
    dphidrc(ntheta) = dphidrc(1)
    dphidthc(1:ntheta-1) = 0.5*(dphidth(1:ntheta-1) + dphidth(2:ntheta))
    dphidthc(ntheta) = dphidth(1)

    phi_neo = phineo(ir_loc,:)

    ! vpadHdE is the derivative of F1/F0 with respect to energy at fixed mu (not xi)
    do ie = 1, nenergy
       do ixi = 1, nxi
          do ig = 1, ntheta
             vpadHdE(ig,ixi,ie,:) = sqrt(energy(ie))*xi(ig,ixi)*dHdE(ig,ixi,ie,:) &
                  + (1.-xi(ig,ixi)**2)*dHdxi(ig,ixi,ie,:)/(2.*sqrt(energy(ie)))
          end do
       end do
    end do

    do is = 1, ns
       do ie = 1, nenergy
          do il = 1, nlambda
             ! should be weighted with backdif for consistency -- MAB
             dHdrc(1:ntheta-1,il,ie,1,is) = 0.5*(dHdr(1:ntheta-1,2*nlambda-il,ie,is) &
                  + dHdr(2:ntheta,2*nlambda-il,ie,is))
             dHdrc(1:ntheta-1,il,ie,2,is) = 0.5*(dHdr(1:ntheta-1,il,ie,is) &
                  + dHdr(2:ntheta,il,ie,is))
             dHdthc(1:ntheta-1,il,ie,1,is) = 0.5*(dHdth(1:ntheta-1,2*nlambda-il,ie,is) &
                  + dHdth(2:ntheta,2*nlambda-il,ie,is))
             dHdthc(1:ntheta-1,il,ie,2,is) = 0.5*(dHdth(1:ntheta-1,il,ie,is) &
                  + dHdth(2:ntheta,il,ie,is))
             dHdEc(1:ntheta-1,il,ie,1,is) = 0.5*(dHdE(1:ntheta-1,2*nlambda-il,ie,is) &
                  + dHdE(2:ntheta,2*nlambda-il,ie,is))
             dHdEc(1:ntheta-1,il,ie,2,is) = 0.5*(dHdE(1:ntheta-1,il,ie,is) &
                  + dHdE(2:ntheta,il,ie,is))
             dHdxic(1:ntheta-1,il,ie,1,is) = 0.5*(dHdxi(1:ntheta-1,2*nlambda-il,ie,is) &
                  + dHdxi(2:ntheta,2*nlambda-il,ie,is))
             dHdxic(1:ntheta-1,il,ie,2,is) = 0.5*(dHdxi(1:ntheta-1,il,ie,is) &
                  + dHdxi(2:ntheta,il,ie,is))
             vpadHdEc(1:ntheta-1,il,ie,1,is) = 0.5*(vpadHdE(1:ntheta-1,2*nlambda-il,ie,is) &
                  + vpadHdE(2:ntheta,2*nlambda-il,ie,is))
             vpadHdEc(1:ntheta-1,il,ie,2,is) = 0.5*(vpadHdE(1:ntheta-1,il,ie,is) &
                  + vpadHdE(2:ntheta,il,ie,is))
             hneoc(1:ntheta-1,il,ie,1,is) = 0.5*(hneo(ir_loc,1:ntheta-1,2*nlambda-il,ie,is) &
                  + hneo(ir_loc,2:ntheta,2*nlambda-il,ie,is))
             hneoc(1:ntheta-1,il,ie,2,is) = 0.5*(hneo(ir_loc,1:ntheta-1,il,ie,is) &
                  + hneo(ir_loc,2:ntheta,il,ie,is))
          end do
       end do
    end do
    dHdrc(ntheta,:,:,:,:) = 0.0 ; dHdthc(ntheta,:,:,:,:) = 0.0 ; vpadHdEc(ntheta,:,:,:,:) = 0.0
    dHdEc(ntheta,:,:,:,:) = 0.0 ; dHdxic(ntheta,:,:,:,:) = 0.0 ; hneoc(ntheta,:,:,:,:) = 0.0
    dphidrc(ntheta) = 0.0 ; dphidthc(ntheta) = 0.0

    deallocate (xi, emax, chebyp1, chebyp2, legp, dHdr, dHdth, dHdxi, dHdE, vpadHdE, hneo, forbid)
    deallocate (coefs, drad_neo, dens_neo, temp_neo)
    deallocate (dphidr, dl_over_b, transport)
    deallocate (pflx, qflx, vflx, upar1, qparB)
    deallocate (upar, qpar, uparB)
    deallocate (upar2, qpar2)
    deallocate (afac, bfac)

  end subroutine get_lowflow_terms

  function zfnc (enrgy, enrgymax)

    implicit none

    real, intent(in) :: enrgy, enrgymax
    real :: zfnc

    zfnc = 2.*sqrt(enrgy/enrgymax)-1.

    return

  end function zfnc

  ! knd = 1 (2) for cheb polys of first (second) kind
  subroutine chebyshev (x, chebyp, knd)

    implicit none

    integer, intent (in) :: knd
    real, intent (in) :: x
    real, dimension (0:), intent (out) :: chebyp

    integer :: n, idx

    n = size(chebyp)-1

    chebyp(0) = 1.0

    if (knd == 2) then
       chebyp(1) = 2.*x
    else
       chebyp(1) = x
    end if
    do idx = 2, n
       chebyp(idx) = 2.*x*chebyp(idx-1) - chebyp(idx-2)
    end do

  end subroutine chebyshev

  subroutine legendre (x, legp)

    implicit none

    real, intent (in) :: x
    real, dimension (0:), intent (out) :: legp

    integer :: n, idx

    n = size(legp)-1

    legp(0) = 1.0
    legp(1) = x

    do idx = 2, n
       legp(idx) = ((2.*idx-1.)*x*legp(idx-1) + (1.-idx)*legp(idx-2))/idx
    end do

  end subroutine legendre

  subroutine get_H (gjk, legdre, chebyshv, h, phi, ig, is)

    use species, only: spec

    implicit none

    real, dimension (0:,0:), intent (in) :: gjk
    real, dimension (:,0:), intent (in) :: legdre, chebyshv
    real, intent (in) :: phi
    real, dimension (:,:), intent (out) :: h
    integer, intent (in) :: ig, is

    integer :: ix, ij, ik

    h = 0.0

    do ix = 1, size(h,1)
       do ik = 0, size(gjk,2)-1
          do ij = 0, size(gjk,1)-1
             h(ix,:) = h(ix,:) + chebyshv(:,ik)*legdre(ix,ij)*gjk(ij,ik)
          end do
       end do
    end do

    ! F_1s = H_1s - Z_s*e*Phi_1/T_s
    ! want Z_s*e*Phi/T_s, but phi from neo is e*Phi/Tnorm
    ! at center (GS2) radius, this is a simple multiply by Z_s * Tref/Ts
    ! note that this assumes the Tref used in NEO is the same as that used in GS2
    h = h - spec(is)%zt*phi

  end subroutine get_H

  subroutine get_dHdxi (gjk, legdre, chebyshv, x, dH)

    implicit none

    real, dimension (0:,0:), intent (in) :: gjk
    real, dimension (:,0:), intent (in) :: legdre
    real, dimension (:), intent (in) :: x
    real, dimension (:,0:), intent (in) :: chebyshv
    real, dimension (:,:), intent (out) :: dH

    integer :: ij, ik, ix

    dH = 0.0

    ! calculate dH/dxi = sum_{j,k} g_{j,k} * T_{k}(z) * d(P_{j}(xi))/d(xi), where z=2*sqrt(E/EMAX)-1
    ! note that dP_{j}/dxi = (j+1)*(P_{j+1} - xi*P_{j})/(xi^2-1)
    do ix = 1, size(x)
       do ik = 0, size(gjk,2)-1
          do ij = 0, size(gjk,1)-1
             dH(ix,:) = dH(ix,:) + (chebyshv(:,ik)*(ij+1)*(legdre(ix,ij+1)-x(ix)*legdre(ix,ij))*gjk(ij,ik))/(x(ix)**2-1.)
          end do
       end do
    end do

  end subroutine get_dHdxi

  subroutine get_dHdE (gjk, legdre, chebyshv1, chebyshv2, x, xmax, dH)

    implicit none

    real, dimension (0:,0:), intent (in) :: gjk
    real, dimension (:,0:), intent (in) :: legdre
    real, dimension (:), intent (in) :: x
    real, intent (in) :: xmax
    real, dimension (:,0:), intent (in) :: chebyshv1, chebyshv2
    real, dimension (:,:), intent (out) :: dH

    integer :: ij, ik, ix
    real :: cfac

    dH = 0.0

    ! dH/dE at fixed xi is sum_{j,k} c_{j,k}*P_j(xi)*d(T_{n}(z))/dz / sqrt(E*EMAX)
    ! z=2*sqrt(E/EMAX)-1
    ! note that d(T_{k})/dz = k*U_{k-1}, where U_k is the Chebyshev polynomial of the 2nd kind
    do ix = 1, size(x)
       do ik = 0, size(gjk,2)-1
          if (ik==0) then
! cfac nonzero here if hneo = F_1 instead of F_1/F_0
!             cfac = -chebyshv1(ix,ik)
             cfac = 0.0
          else
! extra chebyshv1 term not needed since hneo = F_1/F_0 instead of F_1
!             cfac = ik/sqrt(xmax*x(ix))*chebyshv2(ix,ik-1)-chebyshv1(ix,ik)
             cfac = ik/sqrt(xmax*x(ix))*chebyshv2(ix,ik-1)
          end if
          do ij = 0, size(gjk,1)-1
             dH(:,ix) = dH(:,ix) + (legdre(:,ij)*cfac*gjk(ij,ik))
          end do
       end do
    end do

  end subroutine get_dHdE

  ! returns radial derivative at center of 3 points
  subroutine get_radgrad (h, rad, ir, dh)

    implicit none

    real, dimension (:), intent (in) :: h, rad
    integer, intent (in) :: ir
    real, intent (out) :: dh

    dh = (h(ir+1)-h(ir-1))/(rad(ir+1)-rad(ir-1))

  end subroutine get_radgrad

  subroutine get_thgrad (h, th, dh)

    implicit none

    real, dimension (:), intent (in) :: h, th
    real, dimension (:), intent (out) :: dh

    integer :: ig, nth

    nth = size(th)

    do ig = 2, nth-1
       dh(ig) = (h(ig+1)-h(ig-1))/(th(ig+1)-th(ig-1))
    end do
    ! note that H_neo is periodic in theta
    dh(1) = (h(2)-h(nth))/(2.*(th(2)-th(1)))
    dh(nth) = (h(1)-h(nth-1))/(2.*(th(nth)-th(nth-1)))

  end subroutine get_thgrad

  subroutine get_gradH (dgjk, dphi, legdre, chebyshv, dH)

    implicit none

    real, dimension (0:,0:), intent (in) :: dgjk
    real, intent (in) :: dphi
    real, dimension (:,0:), intent (in) :: legdre
    real, dimension (:,0:), intent (in) :: chebyshv
    real, dimension (:,:), intent (out) :: dH

    integer :: ij, ik, ix

    dH = 0.0

    ! calculate dH/drho or dH/dtheta at fixed (xi,E) as sum_{j,k} d(g_{j,k})/drho * T_{k}(z) * P_{j}(xi)
    ! similarly for dH/dtheta
    ! where z=2*sqrt(E/EMAX)-1
    do ix = 1, size(legdre,1)
       do ik = 0, size(dgjk,2)-1
          do ij = 0, size(dgjk,1)-1
             dH(ix,:) = dH(ix,:) + chebyshv(:,ik)*legdre(ix,ij)*dgjk(ij,ik)
          end do
       end do
    end do

    dH = dH - dphi

  end subroutine get_gradH

  subroutine read_neocoefs (theta, nspec_neo, nenergy_neo, nxi_neo, nrad_neo, ir_neo)

    use mp, only: proc0, broadcast
    use splines, only: lf_spline
    use file_utils, only: get_unused_unit
    use constants, only: pi,twopi

    implicit none

    real, dimension (:), intent (in) :: theta
    integer, intent (out) :: nspec_neo, nenergy_neo, nxi_neo, nrad_neo, ir_neo

    integer :: is, ik, ij, ig, ir, idx, ntheta, ntheta_neo
    integer :: ptheta_neo, ip
    integer :: ngrid_per_rad
    integer, save :: neo_unit, neof_unit, neophi_unit

    real, dimension (:), allocatable :: tmp, neo_coefs, dum, neo_phi, dneo_phi
    real, dimension (:), allocatable :: theta_neo_ext,neo_coefs_ext, neo_phi_ext, dneo_phi_ext !JPL
    ntheta = size(theta)

    ! for now, set ir_neo by hand, but best to derive it from neo output in future
    ir_neo = 2

    if (proc0) then
       call get_unused_unit (neo_unit)

       ! read in number of grid points from neo's grid.out file
       open (unit=neo_unit, file='neo_grid.out', status="old", action="read")
       read (neo_unit,*) nspec_neo
       read (neo_unit,*) nenergy_neo
       read (neo_unit,*) nxi_neo
       read (neo_unit,*) ntheta_neo
       if (.not. allocated(theta_neo)) allocate (theta_neo(ntheta_neo))
       do ig = 1, ntheta_neo
          read (neo_unit,*) theta_neo (ig)
       end do
       read (neo_unit,*) nrad_neo
       if (.not. allocated(rad_neo)) allocate (rad_neo(nrad_neo))
       if (.not. allocated(dens_neo)) allocate (dens_neo(nrad_neo,nspec_neo))
       if (.not. allocated(temp_neo)) allocate (temp_neo(nrad_neo,nspec_neo))
       do ir = 1, nrad_neo
          read (neo_unit,*) rad_neo(ir), dens_neo(ir,:), temp_neo(ir,:)
       end do
       close (neo_unit)
       dens_neo = dens_neo / spread(dens_neo(ir_neo,:),1,nrad_neo)
       temp_neo = temp_neo / spread(temp_neo(ir_neo,:),1,nrad_neo)
    end if

    call broadcast (nspec_neo)
    call broadcast (nenergy_neo)
    call broadcast (nxi_neo)
    call broadcast (ntheta_neo)
    call broadcast (nrad_neo)
    if (.not. allocated(theta_neo)) allocate (theta_neo(ntheta_neo))
    if (.not. allocated(rad_neo)) allocate (rad_neo(nrad_neo))
    if (.not. allocated(dens_neo)) allocate (dens_neo(nrad_neo,nspec_neo))
    if (.not. allocated(temp_neo)) allocate (temp_neo(nrad_neo,nspec_neo))
    call broadcast (theta_neo)
    call broadcast (rad_neo)
    do ir = 1, nrad_neo
       call broadcast (dens_neo)
       call broadcast (temp_neo)
    end do

    ! define number of grid points per radial location to make file i/o easier
    ngrid_per_rad = ntheta_neo*(nxi_neo+1)*nenergy_neo*nspec_neo

    if (.not. allocated(tmp)) allocate (tmp(ngrid_per_rad*nrad_neo))
    if (.not. allocated(neo_coefs)) allocate (neo_coefs(ntheta_neo), neo_phi(ntheta_neo), dneo_phi(ntheta_neo))
    if (.not. allocated(dum)) allocate (dum(ntheta))
    if (.not. allocated(coefs)) allocate (coefs(nrad_neo,ntheta,0:nxi_neo,0:nenergy_neo-1,nspec_neo))
    if (.not. allocated(phineo)) allocate (phineo(nrad_neo,ntheta))
    if (.not. allocated(dphidth)) allocate (dphidth(ntheta))
    if (.not. allocated(dcoefsdr)) then
       allocate (dcoefsdr(0:nxi_neo,0:nenergy_neo-1)) ; dcoefsdr = 0.
       allocate (dcoefsdth(ntheta,0:nxi_neo,0:nenergy_neo-1)) ; dcoefsdth = 0.
    end if

    if (proc0) then

       !JPL: if the range of theta grid in GS2 is beyond [-pi:pi](e.g. "nperiod > 1"),
       !     extend the theta grid in NEO ([-pi:pi]) to the periodic one in the theta range of GS2.
       ptheta_neo=ceiling((maxval(theta)+pi)/twopi)  !the period of 2pi in theta range
       if (ptheta_neo .gt. 1) then
          if (.not. allocated(theta_neo_ext)) allocate (theta_neo_ext(ntheta_neo*(2*ptheta_neo-1)))
          if (.not. allocated(neo_coefs_ext)) allocate (neo_coefs_ext(ntheta_neo*(2*ptheta_neo-1)))
          if (.not. allocated(neo_phi_ext)) allocate (neo_phi_ext(ntheta_neo*(2*ptheta_neo-1)))
          if (.not. allocated(dneo_phi_ext)) allocate (dneo_phi_ext(ntheta_neo*(2*ptheta_neo-1)))
       end if

       ! read in H1^{nc} (adiabatic piece of F1^{nc}) from neo's f.out file
       call get_unused_unit (neof_unit)
       open (unit=neof_unit, file='neo_f.out', status="old", action="read")

       do ir = 1, nrad_neo
          read (neof_unit,*) tmp(ngrid_per_rad*(ir-1)+1:ngrid_per_rad*ir)
       end do

       idx = 1
       do ir = 1, nrad_neo
          do is = 1, nspec_neo
             do ik = 0, nenergy_neo-1
                do ij = 0, nxi_neo
                   do ig = 1, ntheta_neo
                      neo_coefs(ig) = tmp(idx)
                      idx = idx+1
                   end do

                   if ( ptheta_neo .gt. 1) then
                      do ip = 1, ntheta_neo*(2*ptheta_neo-1)
                         theta_neo_ext(ip)=theta_neo(mod((ip-1),ntheta_neo)+1)+twopi*(int((ip-1)/ntheta_neo)-(ptheta_neo-1))
                         neo_coefs_ext(ip)=neo_coefs(mod((ip-1),ntheta_neo)+1)
                      end do
                      call lf_spline (theta_neo_ext, neo_coefs_ext, theta, coefs(ir,:,ij,ik,is), dum)  !JPL
                   else
                   ! need to interpolate coefficients from neo's theta grid to gs2's
                      call lf_spline (theta_neo, neo_coefs, theta, coefs(ir,:,ij,ik,is), dum)
                   end if
                end do
             end do
          end do
       end do

       close (neof_unit)
    end if

    deallocate (tmp)

    if (.not. allocated(tmp)) allocate (tmp(ntheta_neo*nrad_neo))

    if (proc0) then
       ! read in phi1^{nc} from neo's phi.out file
       call get_unused_unit (neophi_unit)
       open (unit=neophi_unit, file='neo_phi.out', status="old", action="read")

       do ir = 1, nrad_neo
          read (neophi_unit,*) tmp((ir-1)*ntheta_neo+1:ir*ntheta_neo)
       end do

       idx = 1
       do ir = 1, nrad_neo
          do ig = 1, ntheta_neo
             neo_phi(ig) = tmp(idx)
             idx = idx+1
          end do

          ! TMP FOR TESTING -- MAB
!          neo_phi = 0.

          ! need to interpolate coefficients from neo's theta grid to gs2's
          if ( ptheta_neo .gt. 1) then
             do ip = 1, ntheta_neo*(2*ptheta_neo-1)
                neo_phi_ext(ip)=neo_phi(mod((ip-1),ntheta_neo)+1)
             end do
             call lf_spline (theta_neo_ext, neo_phi_ext, theta, phineo(ir,:), dum)  !JPL
          else
             call lf_spline (theta_neo, neo_phi, theta, phineo(ir,:), dum)
          end if

          ! at central radius, calculate dphi/dth and interpolate onto gs2 grid
          if (ir == 2) then
             if ( ptheta_neo .gt. 1) then
                do ip = 1, ntheta_neo*(2*ptheta_neo-1)
                   dneo_phi_ext(ip)=dneo_phi(mod((ip-1),ntheta_neo)+1)
                end do
                call get_thgrad (neo_phi_ext, theta_neo_ext, dneo_phi_ext)
                call lf_spline (theta_neo_ext, dneo_phi_ext, theta, dphidth, dum)
             else
                call get_thgrad (neo_phi, theta_neo, dneo_phi)
                call lf_spline (theta_neo, dneo_phi, theta, dphidth, dum)
             end if
          end if
       end do

       close (neophi_unit)

    end if

    call broadcast (dphidth)
    do ir = 1, nrad_neo
       call broadcast (phineo(ir,:))
    end do
    do ir = 1, nrad_neo
       do is = 1, nspec_neo
          do ik = 0, nenergy_neo-1
             do ij = 0, nxi_neo
                call broadcast (coefs(ir,:,ij,ik,is))
             end do
          end do
       end do
    end do

    deallocate (tmp, neo_coefs, theta_neo, neo_phi, dneo_phi, dum)
    if (proc0 .and. (ptheta_neo .gt. 1)) then
       deallocate (neo_coefs_ext, theta_neo_ext, neo_phi_ext, dneo_phi_ext)
    endif
  end subroutine read_neocoefs

end module lowflow


text:
	SUBROUTINE DSTANM  ( iret )
C************************************************************************
C* DSTANM								*
C* 									*
C* This subroutine starts a new animation sequence.			*
C* 									*
C* DSTANM  ( IRET )							*
C*									*
C* Output parameters:							*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* A. Chang/EAI	 	12/93						*
C* S. Jacobs/NMC         3/94           Renamed DSPIXM to DSTANM        *
C* S. Jacobs/NMC	 6/94		Added call to HSTANM		*
C************************************************************************
	INCLUDE		'ERROR.PRM'
C------------------------------------------------------------------------
	iret = NORMAL
C
C*	Call driver to start the sequence.
C
	CALL HSTANM ( iret )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      program lagrit_main
C
C #####################################################################
C
C     PURPOSE -
C
C        start up lagrit
C
C     INPUT ARGUMENTS -
C
C        None
C
C     OUTPUT ARGUMENTS -
C
C        None
C
C     CHANGE HISTORY -
C
C        $Log: lagrit_main.f,v $
C        Revision 2.00  2007/11/05 19:48:30  spchu
C        Import to CVS and change file name from adrivgen.f to lagrit_main.f
C
CPVCS       Rev 1.13   31 Jan 2019 13:29:00   drl
CPVCS    added basic command line parsing
CPVCS
CPVCS       Rev 1.12   09 Mar 2000 10:07:06   dcg
CPVCS    changes for no storage block version of lagrit
CPVCS
CPVCS       Rev 1.11   Fri Oct 22 11:00:36 1999   dcg
CPVCS    get rid of subroutine fiximts
CPVCS
CPVCS       Rev 1.10   Mon Feb 22 16:16:10 1999   dcg
CPVCS    rewrite of command processing to allow for recursion
C ######################################################################
C
      implicit none
      integer ierror_return, i
      character(len=32) :: arg,arg_value,mode,logfile,outfile
C
      mode = 'noisy'
      logfile = ' '
      outfile = ' '

      i = 0
      do
C        Get command line flag
         call get_command_argument(i,arg)
         if (len_trim(arg) == 0) exit

C        Get flag value - will return empty string if it doesn't exist
         call get_command_argument(i+1,arg_value)

C        Verify that flag value is not empty and is not another flag
         if ((arg_value(1:1).ne.'-').and.(arg_value(1:1).ne.' ')) then
            i = i + 1

            select case (trim(arg))
            case ('-log')
               logfile = trim(arg_value)
            case ('-out')
               outfile = trim(arg_value)
            case ('-mode')
               mode = trim(arg_value)
               if (mode(1:6).ne.'silent') then
                  if (mode(1:5).ne.'noisy') then
                     print*,'ERROR: Invalid option for -mode: ',mode
                     stop
                  endif
               endif
            end select
         endif

         i = i + 1
      end do

      ierror_return = 0
      call initlagrit(trim(mode),logfile,outfile)
      call control_command_lg(ierror_return)

      stop
      end program lagrit_main

max_stars_repo_name: daniellivingston/LaGriT-Windows

max_issues_repo_name: daniellivingston/LaGriT

max_forks_repo_name: daniellivingston/LaGriT


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 09.02.10/04                           *
C  *    TEST TITLE : Error indicator = 2004                *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /ERRINF/ ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR,
     1                USRERR, ERRSAV,      FUNSAV,      FILSAV,
     2                EFCNT, EFID
      INTEGER         ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR(10),
     1                USRERR, ERRSAV(200), FUNSAV(200), FILSAV(200),
     2                EFCNT, EFID(100)
      COMMON /ERRCHR/ CURCON,     ERRSRS,    ERRMRK,    ERFLNM,
     1                CONTAB
      CHARACTER       CURCON*200, ERRSRS*40, ERRMRK*20, ERFLNM*80,
     1                CONTAB(40)*150

      INTEGER     IA210(2,10), SPATH(2,1), EIS(2), EES(1),IDUM1,ISIZ
      INTEGER     I1(2),I2(2),I3(2),I4(2),I5(2),I6(2),I7(2),I8(2)
      REAL        PX(1),PY(1), RA(3)
      DATA        PX/1./,PY/0.5/

      CALL INITGL ('09.02.10/04')

      CALL XPOPPH (ERRFIL, MEMUN)

      CURCON = 'input parameter size is out of range'
      CALL SETVS ('2004', EXPERR,EXPSIZ)
      ERRSRS = '4'

      CALL POPST (100)
      CALL PTX (0.5,0.5,'TEXT100')
      CALL PEXST (101)
      CALL PCLST
      CALL POPST (101)
      CALL PPM (1, PX, PY)
      CALL PCLST

      CALL SETVS ('1,3', EIS, ISIZ)
      CALL SETVS ('2', EES, ISIZ)
      CALL RELS (100,1,1,-1,EIS,1,EES,ERRIND,IDUM1,IDUM1)
      CALL RELS (100,1,1, 2,EIS,-1,EES,ERRIND,IDUM1,IDUM1)

      CALL SETVS ('100,1', SPATH, ISIZ)
      CALL RISS3 (0.,0.1,0.,0.1,-1, SPATH,0,0,
     1            0,I1,I2,I3,I4, 0, I5,I6,I7,I8,
     2            10, ERRIND, IDUM1, IA210)

      CALL POPWK (WKID, CONID, WTYPE)
      CALL RQCR  (WKID, 1, -1, 0, ERRIND, IDUM1, RA)
      CALL PCLWK (WKID)

      CALL ENDIT

      END


text:
      subroutine absswb(adhfs,adhfsx,advfs,advfsx,avcnst,axhfs,
     $ axhfsx,axlks,axlksx,axvfs,axvfsx,cdrs,cdrtw,cdrsx,cdrw,
     $ csts,dhfs,dvfs,eps100,ibswx,iindx1,iodb,ipch,ipchmx,ipcv,
     $ ipcvmx,jcsort,jflag,jsflag,kbt,kmax,mosp,mtb,narn1,narn2,
     $ narxmx,narxt,nbasp,nbaspd,nbaspx,nbt,nbtmax,nbw,ndrs,ndrsmx,
     $ ndrsr,ndrsrx,ndrsx,nelect,nhydr,nodbmx,no2gaq,noutpt,nst,
     $ nstmax,nsts,nstsmx,nstsr,nswtch,ntpr,ntprmx,nttyo,presg,
     $ press,qbassw,qbswx,tempc,uspec,uzvec1,weight,xhfs,xvfs,xlks)
c
c     This subroutine carries out automatic basis switching to optimize
c     optimize the basis set after Newton-Raphson iteration has
c     converged at the current point of reaction progress. Such
c     optimization tends to improve the code numerics (finite
c     difference stability, matrix conditioning, etc.) This subroutine
c     is similar to EQLIB/absswa.f, which carries out automatic basis
c     switching as a pre-Newton-Raphson optimization technique.
c
c     This subroutine is called by:
c
c       EQ6/path.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c
c     Principal output:
c
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer ipchmx,ipcvmx,kmax,narxmx,nbtmax,ndrsmx,nodbmx,nstmax,
     $ nstsmx,ntprmx
c
      integer noutpt,nttyo
c
      integer ibswx(nbtmax),iindx1(kmax),iodb(nodbmx),jcsort(nstmax),
     $ jflag(nstmax),jsflag(nstmax),narxt(ntprmx),nbasp(nbtmax),
     $ nbaspd(nbtmax),nbaspx(nbtmax),ndrs(ndrsmx),ndrsr(2,nstmax),
     $ ndrsrx(2,nstmax),ndrsx(ndrsmx),nsts(nstsmx),nstsr(2,nstmax)
c
      integer ipch,ipcv,kbt,narn1,narn2,nbt,nbw,nelect,nhydr,no2gaq,
     $ nst,nswtch,ntpr
c
      logical qbassw,qbswx
c
      character*48 uspec(nstmax),uzvec1(kmax)
c
      real*8 adhfs(narxmx,ntprmx,ipchmx,nstmax),
     $ adhfsx(narxmx,ntprmx,ipchmx,nstmax),
     $ advfs(narxmx,ntprmx,ipcvmx,nstmax),
     $ advfsx(narxmx,ntprmx,ipcvmx,nstmax),
     $ axhfs(narxmx,ntprmx,nstmax),axhfsx(narxmx,ntprmx,nstmax),
     $ axlks(narxmx,ntprmx,nstmax),axlksx(narxmx,ntprmx,nstmax),
     $ axvfs(narxmx,ntprmx,nstmax),axvfsx(narxmx,ntprmx,nstmax),
     $ cdrs(ndrsmx),cdrsx(ndrsmx),cdrtw(nstmax),cdrw(nstmax),
     $ csts(nstsmx),dhfs(ipchmx,nstmax),dvfs(ipcvmx,nstmax),
     $ mosp(nstmax),mtb(nbtmax),weight(nstmax),xhfs(nstmax),
     $ xlks(nstmax),xvfs(nstmax)
c
      real*8 avcnst,presg,press,tempc
c
      real*8  eps100
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer jlen1,jlen2,jlen3,krow,nb,nb1,nb2,ns,nse,nsi,nsj,ns1,ns2
c
      character*56 usp156,usp256,usp356
c
      real*8 api,apj,cx,fx1,fx2,rx,wb1,wb2,wsi
c
      real*8 coefdr,coefst
c
c-----------------------------------------------------------------------
c
c     Initialize the number of basis switches made.
c
      nswtch = 0
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Pick candidates for automatic basis switching. The indices of
c     candidates are stored in the array ibswx. Note that the method
c     employed here is similar to that in EQLIB/abswpk.f. However,
c     there are some differences. Here, basis switching is not being
c     used as a method of pre-Newton-Raphson optimization. Hence, there
c     are no residual functions used in defining criteria for making
c     a switch. the relative contributions to mass balances are used
c     instead.
c
      qbswx = .false.
      do nb = 1,nbt
        ibswx(nb) = 0
        nse = nbaspd(nb)
        nsj = nbasp(nb)
        if (nse.ne.narn1 .and. nse.ne.nhydr .and. nse.ne.no2gaq
     $    .and.nse.ne.nelect) then
c
c         If the current species is not H2O, H+, O2(g,aq), or e-,
c         consider a switch.
c
c         Get weights for the current mass balance.
c
          do ns = narn1,narn2
            weight(ns) = coefst(csts,nsts,nstsmx,nstsr,nb,ns,nstmax)
          enddo
c
c         Screen out species in the mass balance that are not linked
c         by the current reaction with the current basis species.
c         Do this by seting the corresponding weights to zero.
c
          do ns = narn1,narn2
            if (jflag(ns) .eq. 30) then
c
c             Calling sequence substitutions:
c               nsj for nse
c
              cx = coefdr(cdrs,ndrs,ndrsmx,ndrsr,nsj,ns,nstmax)
              if (cx .eq. 0.) weight(ns) = 0.
            endif
          enddo
c
c         Find a candidate.
c
          call fbassw(jcsort,jflag,mosp,narn1,narn2,nse,nsi,nsj,
     $    nstmax,weight,wsi)
c
c         If a candidate was found, apply a filter.
c
          if (nsi .gt. 0) then
c
c           Use the current candidate species if its contribution
c           to the mass balance exceeds that of the current basis
c           species by a reasonable factor, here 10.
c
            api = weight(nsi)*mosp(nsi)
            apj = weight(nsj)*mosp(nsj)
            rx = api/apj
            if (rx .gt. 10.) then
              ibswx(nb) = nsi
              qbswx = .true.
            endif
          endif
        endif
      enddo
c
c     If there are no candidates for basis switching, exit the
c     present subroutine.
c
      if (.not.qbswx) go to 999
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (iodb(1) .ge. 2) then
        write (noutpt,1000)
 1000   format(/5x,'--- Candidate Basis Switches ---',/)
        do nb = 1,nbt
          nse = nbaspd(nb)
          nsj = nbasp(nb)
          nsi = ibswx(nb)
          if (nsi .gt. 0) then
c
c           Calling sequence substitutions:
c             jlen1 for jlen
c             uspec(nsj) for unam48
c             usp156 for uspn56
c
            call fmspnx(jlen1,uspec(nsj),usp156)
c
c           Calling sequence substitutions:
c             jlen2 for jlen
c             uspec(nsi) for unam48
c             usp256 for uspn56
c
            call fmspnx(jlen2,uspec(nsi),usp256)
c
c           Calling sequence substitutions:
c             jlen3 for jlen
c             uspec(nse) for unam48
c             usp356 for uspn56
c
            call fmspnx(jlen3,uspec(nse),usp356)
c
            write (noutpt,1010) usp156(1:jlen1),usp256(1:jlen2),
     $      usp356(1:jlen3)
 1010       format(/3x,'Could replace ',a,' in the active basis set',
     $      ' with',/3x,a,' as the species associated with the mass',
     $      ' balance',/3x,'of ',a,'.')
          endif
        enddo
        write (noutpt,1020)
 1020   format(1x)
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Resolve any conflicts in candidate basis switches. The algorithm
c     used here is similar to but not identical to that employed by
c     EQLIB/gabswx.f, which resolves conflicts in the context of using
c     automatic basis switching as a form of pre-Newton-Raphson
c     optimization.
c
      do nb1 = 1,nbt - 1
        ns1 = ibswx(nb1)
        if (ns1 .gt. 0) then
          do nb2 = nb1 + 1,nbt
            ns2 = ibswx(nb2)
            if (ns2 .eq. ns1) then
c
c             Calling sequence substitutions:
c               nb1 for nb
c
              wb1 = coefst(csts,nsts,nstsmx,nstsr,nb1,ns,nstmax)
              fx1 = wb1*mosp(ns1)/mtb(nb1)
c
c             Calling sequence substitutions:
c               nb2 for nb
c
              wb2 = coefst(csts,nsts,nstsmx,nstsr,nb2,ns,nstmax)
              fx2 = wb2*mosp(ns2)/mtb(nb2)
              if (fx1 .gt. fx2) then
                ibswx(nb2) = 0
              else
                ibswx(nb1) = 0
                go to 100
              endif
            endif
          enddo
        endif
  100   continue
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      nswtch = 0
      do nb = 1,nbt
        nsi = ibswx(nb)
        if (nsi .gt. 0) nswtch = nswtch + 1
      enddo
c
      if (nswtch .le. 0) then
        if (iodb(1) .ge. 2) write (noutpt,1070)
 1070   format(/10x,'No switches will be made.',/)
        go to 999
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (iodb(1) .ge. 2) then
        write (noutpt,1100)
 1100   format(/5x,'--- Automatic Basis Switches ---',/)
        do nb = 1,nbt
          nse = nbaspd(nb)
          nsj = nbasp(nb)
          nsi = ibswx(nb)
          if (nsi .gt. 0) then
c
c           Calling sequence substitutions:
c             jlen1 for jlen
c             uspec(nsj) for unam48
c             usp156 for uspn56
c
            call fmspnx(jlen1,uspec(nsj),usp156)
c
c           Calling sequence substitutions:
c             jlen2 for jlen
c             uspec(nsi) for unam48
c             usp256 for uspn56
c
            call fmspnx(jlen2,uspec(nsi),usp256)
c
c           Calling sequence substitutions:
c             jlen3 for jlen
c             uspec(nse) for unam48
c             usp356 for uspn56
c
            call fmspnx(jlen3,uspec(nse),usp356)
c
            write (noutpt,1110) usp156(1:jlen1),usp256(1:jlen2),
     $      usp356(1:jlen3)
 1110       format(/3x,'Will replace ',a,' in the active basis set',
     $      ' with',/3x,a,' as the species associated with the mass',
     $      ' balance',/3x,'of ',a,'.')
          endif
        enddo
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      call autosw(adhfs,adhfsx,advfs,advfsx,axhfs,axhfsx,
     $ axlks,axlksx,axvfs,axvfsx,cdrs,cdrsx,eps100,ibswx,iindx1,
     $ ipch,ipchmx,ipcv,ipcvmx,jflag,jsflag,kbt,kmax,narn1,narxmx,
     $ nbasp,nbaspd,nbaspx,nbt,nbtmax,nbw,ndrs,ndrsmx,ndrsx,ndrsr,
     $ ndrsrx,noutpt,nst,nstmax,ntprmx,nttyo,qbassw,uspec,uzvec1)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Recompute the cdrw array.
c
      call gcdrw(cdrs,cdrw,narn1,ndrs,ndrsmx,ndrsr,nst,nstmax)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Recompute the cdrtw array.
c
      call gcdrtw(cdrs,cdrtw,narn1,narn2,ndrs,ndrsmx,ndrsr,
     $ nelect,no2gaq,nst,nstmax)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Update the thermodynamic data to correspond to the new
c     active basis set. First, recompute the log K, etc., data for
c     the various reactions.
c
      call evdatr(adhfs,advfs,axhfs,axlks,axvfs,dhfs,dvfs,
     $ ipch,ipchmx,ipcv,ipcvmx,narxmx,narxt,nst,nstmax,ntpr,ntprmx,
     $ tempc,xhfs,xlks,xvfs)
c
c     Then make pressure corrections to these thermodynamic data.
c
      if (ipcv .ge. 0) then
        call pcorrx(avcnst,dhfs,dvfs,ipch,ipchmx,ipcv,ipcvmx,
     $  nbasp,nbt,nbtmax,ndrsr,nst,nstmax,presg,press,xhfs,
     $  xlks,xvfs)
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
  999 continue
c
      end


text:
!===============================================================================
!===============================================================================
MODULE os_mod
  USE io_unit_mod
  USE omp_mod
  USE mpi_mod
  USE omp_timer_mod
  implicit none
  private
  type, public:: os_t
  contains
    procedure, nopass:: mkdir
    procedure, nopass:: mkdir_no
  end type
  type(os_t), public:: os
  real, parameter:: ms=100.
CONTAINS

SUBROUTINE mkdir (dir)
  character(len=*):: dir
  real(8):: start
  logical:: exists
  !.............................................................................
  if (mpi%rank==0) then
    call mkdir_real (dir)
  else
    start = wallclock()
    do while (wallclock()-start < 10.0)
#ifdef __INTEL_COMPILER
      inquire (directory=dir, exist=exists)
#else
      inquire (file=dir, exist=exists)
#endif
      if (exists) exit
      call mpi%delay (ms=ms)
    end do
    if (.not.exists) then
      write (io_unit%output,*) 'WARNING: rank',mpi%rank,' cannot see',dir
    end if
  end if
END SUBROUTINE mkdir

!===============================================================================
!> Utility to create a new directory
!===============================================================================
SUBROUTINE mkdir_real (dir)
#ifdef __INTEL_COMPILER
  USE ifport
#else
  use iso_c_binding
  interface
    function fmkdir(path,mode) bind(c,name="mkdir")
      use iso_c_binding
      integer(c_int) :: fmkdir
      character(kind=c_char,len=1) :: path(*)
      integer(c_int16_t), value :: mode
    end function fmkdir
  end interface
  integer :: i, iter
#endif
  character(len=*):: dir
  logical:: exists
  !.............................................................................
#ifdef __INTEL_COMPILER
  inquire (directory=dir, exist=exists)
  if (exists) return
  if (.not.io_unit%do_validate) then
    exists = MakeDirQQ (dir)
    if (exists) then
      write(stderr,*) ' Intel created directory '//dir
    else
      write(stderr,'(2(a,i5,2x),a)') 'rank:', mpi%rank, 'thread:', omp%thread, &
        ' WARNING: Intel failed to create directory '//dir
    end if
  end if
#else
  !$omp critical (system_cr)
  inquire (file=dir, exist=exists)
  if (.not.exists) then
    do iter=1,10
      i = fmkdir(dir, int(o'772',c_int16_t))
      call mpi%delay (ms=ms)
      inquire (file=dir, exist=exists)
      if (exists) then
        if (.not.io_unit%do_validate) &
          write(stderr,*) ' C-binding call created directory '//dir
        exit
      else
        if (.not.io_unit%do_validate) &
          write(stderr,'(2(a,i5,2x),a)') 'rank:', mpi%rank, 'thread:', omp%thread, &
            ' WARNING: C-binding call failed to create directory '//dir
      end if
      call system ('mkdir -p '//dir)
      call mpi%delay (ms=ms)
      inquire (file=dir, exist=exists)
      if (exists) then
        if (.not.io_unit%do_validate) &
          write(stderr,*) ' system call created directory '//dir
        exit
      else
        if (.not.io_unit%do_validate) &
          write(stderr,'(2(a,i5,2x),a)') 'rank:', mpi%rank, 'thread:', omp%thread, &
          ' WARNING: system call failed to create directory '//dir
      end if
    end do
  end if
  !$omp end critical (system_cr)
#endif
END SUBROUTINE mkdir_real

!===============================================================================
!> Make sure the next directory is created ahead of time
!===============================================================================
SUBROUTINE mkdir_no (iout)
  integer:: iout
  character(len=120):: filename
  integer, save:: checked=-1
  !-----------------------------------------------------------------------------
  !$omp critical (mkdir_no_cr)
  if (iout > checked) then
    write (filename,'(a,i5.5,"/")') trim(io_unit%outputname), iout
    if (io_unit%iodir/=iout) call os%mkdir (trim(filename))
    io_unit%iodir = iout
    checked = iout
  end if
  !$omp end critical (mkdir_no_cr)
END SUBROUTINE mkdir_no
END MODULE os_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch