-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathevolvegpd2.f
19583 lines (14356 loc) · 452 KB
/
evolvegpd2.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
C**********************************************************************
C* !!!!!!!!!!!!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!!!!!!! *
C* *
C* If you read this you might have already altered the program in an *
C* unforseen way!!!! Delete it and download it again!!!!!! *
C* *
C* If you want to study the programm, please copy it to another file *
C* name and do not use it for any scientific work!!!!!!!! *
C* *
C* Please contact us, if you want to seriously study the code in more *
C* detail! *
C* *
C**********************************************************************
BLOCK DATA DATPDF
C
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
CHARACTER*10 NAMPDF, NMHDRN, NAMQRK
CHARACTER*12 MRSFLN
LOGICAL LSTX
C
PARAMETER (Z = 1D -10, ZZ = 1D -20)
C MxF is the maximum number of quark flavors
C MxAdF (= 1 gluon + 1 singlet quark + Additional fake flavors
C to save space for storing info on evolution)
PARAMETER (MXX = 1050, MxQ = 25, MxF = 6, MxAdF = 6)
C PARAMETER (MxPN = MxF * 2 + MxAdF)
PARAMETER (MxPN = MxF * 2 + 2)
PARAMETER (MxQX= MxQ * MXX, MxPQX = MxQX * MxPN)
PARAMETER (MxPDF = 20, MXHDRN = 6)
C
COMMON / IOUNIT / NIN, NOUT, NWRT
Common / PdCntrl/ LPrt, LDbg
COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX, DEL,IV
COMMON / QARAY1 / QINI,QMAX, QV(0:MxQ),TV(0:MxQ), NT,JT,NG
COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx
COMMON / PEVLDT / UPD(MXPQX), KF, Nelmt
COMMON / PARPRZ / XMNPRZ(9), QMNPRZ(9), ALF(9),
> NfLPRZ(9), NfAL(9), MORD(9)
COMMON / INITIA / AN(-MxF : MxF), AI(-MxF : MxF),
> BI(-MxF : MxF), CI(-MxF : MxF)
C
COMMON / NAMPDF / NAMPDF(MxPDF), NMHDRN(-1 : MxHDRN),
> NAMQRK(-MxF : MxF)
COMMON / MRSDAT / MRSFLN (3)
C
Data LPrt, LDbg / 1, 0 /
DATA QINI, QMAX, XMIN, XCR / 1.9, 1.001D2, 0.999D-3, 0.3 /
DATA KF, IKNL, IPD0, IHDN / 10, 1, 1, 1 /
DATA NX, NT, JT, LSTX / 40, 6, 1, .FALSE. /
DATA (NfLPRZ(I), I=1,9) / 6,4,5,5,5,6,6,6,6/
DATA (XMNPRZ(I), I=1,9) / 3*1.d-4,6*1.D-5/
DATA (QMNPRZ(I), I=1,9) / 2.25, 2.0,3.2, 2*2.25, 4*2.0 /
C Validity needs to be checked.
DATA (MORD(I), I=1,9) / 2*1, 7*2/
DATA (ALF(I), I=1,9) / 0.2, 0.2, 0.3, 0.19, 0.215,
$ 0.22, 0.225,0.2,0.2/
DATA (NfAL(I), I=1,9) / 3*4, 6*5 /
C
DATA MRSFLN / 'prmz:hmrse', 'prmz:hmrsb', 'GRID3' /
DATA (NAMPDF(I),I=1,11) / 'Ehlq_1 ','Dk-Ow 1 ', 'DFLM_NLLA',
>'KMRS B0 ','MRS92_D0 ','MT90_S1 ','CTEQ1M ', ' ',
>' ', 'QCD_EVL_1','QCD_EVL_2' /
C
DATA (NMHDRN(I), I=-1,5)/ 'AntiProton', 'Neutron', 'Proton',
> 'IsoScalar', 'Pi_Plus', 'Pi_Minus', 'K_Plus' /
C
DATA (NAMQRK(I), I=-6,6) / '-t_Quark', '-b_Quark', '-c_Quark',
> '-s_Quark', '-d_Quark', '-u_Quark', 'Gluon', 'u_Quark',
> 'd_Quark', 's_Quark', 'c_Quark', 'b_Quark', 't_Quark' /
C
C The following data correspond to EHLQ set 1
DATA AI / 7 * Z, 0.5, 0.4, 4 * Z /
DATA BI / 6 * Z, 3.5, 2 * 1.51, 4 * 1./
DATA CI / 3 * 5, 3 * 8.54, 5.9, 3.5, 4.5, 4 * 5./
DATA AN / 3 *ZZ, .081, 2 *.182, 2.62, 1.78, 0.67, 4 *ZZ /
C
C ****************************
END
C
C
FUNCTION PDF (Iset, Ihadron, Iparton, X, Q, DEL, IR)
C ====================================================================
C FUNCTION PDF (Iset, Ihadron, Iparton, X, Q, IR)
C -------------------------------------------------
C Front-end function for parton distributions of pbar, n, p, D, C, Fe
C all tranformed in terms of distribution of proton (PrtPdf)
C Also checks to see if pdf < 0; if so issue warning and set pdf=0.
C Revised 4/4/94 by HLL & WKT:
C PDF retains its name and argument list for compatibility with all
C existing programs which Call this function; PDF switches between
C target hadrons;
C PdfPrt is for proton target; it switches between different Iset's.
C
C Ihadron = -1, 0, 1, 2, 4, 6 :
C pbar, n, p, D, C, Fe
C 5 is "isoscalar-corrected iron" hence = D
C In all cases, adjust the Iparton label
C to convert to the corresponding proton distribution which is
C given in Function PdfPrt;
C ---------------------------------------------------
C If any of the arguments are unphysical, it issues a warning.
C Return error code IR:
C 0 : O.K.
C 1,2: Error from call to PdfPrt (proton parton distribution)
C 3 : X < 0 or X > 1 : unphysical: warning + set PDF = 0;
C 4 : Ihadron out of range;
C 5 : PDF < 0., unphysical: warning + set PDF = 0;
C 6 : Iparton out of range; warning + set PDF = 0;
C -------------------------------------
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
Character Msg*80
PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
COMMON / IOUNIT / NIN, NOUT, NWRT
COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx
DATA
1 HUGE, DUM, D0m, D1p
1 / 1D10, 0.D0, -1D-6, 1.000001D0 /
1 IW1, IW2 / 2*0 /
C
Ier = 0
C ------- Check x-range
IF (X.LE.D0m .or. X.GT.D1p) Then
Ier = 3
Call WARNR(IW1,NWRT,'X out of range in PDF.', 'X', X,
> D0, D1, 1)
TEM = DUM
EndIf
C Check bounds on Ihadron
If (Ihadron.gt.6 .or. Ihadron.lt.-1 .or. Ihadron.eq.3) then
Call WARNI(IW, NWRT,
> 'Only Ihardon=-1,0,1,2,4,5,6 (pbar,n,p,D,*,) are active',
> 'Ihadron', Ihadron, -1,6,1)
Ir = 4
PDF=0.D0
Return
Endif
C Check bounds on Iparton
Jp=Abs(Iparton)
Neff = NFL(Q)
c nfl(q) returns the number of `light' flavors at scale Q - effective
If ( Jp .gt. NEFF) then
C if Jp > Neff, then set PDF=0 and return
Call WARNI(IW, NWRT,
> 'Iparton out of range',
> 'Iparton', Iparton, -Neff, Neff,1)
Ier = 6
PDF = 0D0
Return
Endif
C --- Conversion of Ihadron to proton distributions,
C if necessary ----
If (Jp.eq.1 .or. Jp.eq.2) Then
C For u and d
C Use Isospin symmetry n<->p == u<->d
Ipartner=3-Jp
If (Iparton.lt.0) Ipartner=-Ipartner
If (Ihadron.eq.1) then
C p:
Tem= PdfPrt(Iset, Iparton, X, Q, Ir)
Elseif (Ihadron.eq.-1) then
C pbar:
Tem= PdfPrt(Iset, -Iparton, X, Q, Ir)
Elseif (Ihadron.eq.0) then
C n:
Tem= PdfPrt(Iset, Ipartner, X, Q, Ir)
Elseif (Ihadron.eq.2 .or.Ihadron.eq.4 .or.Ihadron.eq.5) then
C isoscalar: D, C, ...
Tem=( PdfPrt(Iset, Iparton, X, Q, Ir)
> +PdfPrt(Iset, Ipartner, X, Q, Ir) )/2.D0
Elseif (Ihadron.eq.6) then
C Fe:
Tem=( 26.D0* PdfPrt(Iset, Iparton, X, Q, Ir)
> +30.D0* PdfPrt(Iset, Ipartner, X, Q, Ir) )/56.D0
Endif
Else
C For s,c,b,t
Tem= PdfPrt(Iset, Iparton, X, Q, Ir)
Endif
C --- Make sure PDF >= 0 --------
C (unless Iknl<0 - polarized pdf)
c IF (TEM .LT. D0 .and. Iknl .ge. 1 .and. X.GE.DEL) Then
c IF (TEM .LT. D0m .AND. X .LE. 0.9D0) Then
c Call WARNR(IW2,NWRT,'PDF < 0; Set --> 0', 'PDF',TEM,D0,HUGE,1)
c WRITE (NWRT, '(A, 2I5, 2(1PE15.3))')
c > ' ISET, Iparton, X, Q = ', ISET, Iparton, X, Q
c Ier = 5
c EndIf
c TEM = D0
c EndIf
C -------- Return function value and error code
PDF = TEM
IR = Ier
return
C
end
C
FUNCTION PdfPrt (IPDF, LPRTN, XD, QD, IR)
C
C
C This routine gives the parton ( IPARTN ) distribution function inside
C the proton in a chosen evolved or parametrized form ( Ipdf )
C 2. For Ipdf = 10 or 11: steer to results from QCD evolution program;
C 1 - 9: CTEQ parametrizations
C 12 - 20: CTEQ tables or other forms
C > 30 : other parametrizations: see codes below for details.
C It gives the probability distribution, not the momemtum-weighted one.
C
C Return Code: IER = 0: No error
c < 10: PDF fundtion error
C = 10: Ipdf out of range
c > 10: multi-error
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
C
PARAMETER (MXPDF = 20, MXHDRN = 6, MXF = 6, MXPN = MXF*2+2)
C
CHARACTER MSG*75
C
COMMON / IOUNIT / NIN, NOUT, NWRT
C
DATA IW1 / 0 /
IER = 0
Irr = 0
C
Iprtn = LPRTN
x = XD
Q = QD
C ==== Begin of overall Ipdf IfBlock =====
If (Ipdf .eq. 10) Then
C 10 evolve
Tmp = ParDis (Iprtn, X, Q)
Else
Ier = 10
MSG=
>'Ipdf chosen is currently inactive. PdfPrt set equal to zero.'
CALL WARNI (IW1, NWRT, MSG, 'Ipdf', Ipdf, 1, 10, 0)
Tmp = 0.
EndIf
PdfPrt = Tmp
Ir = Ier + Irr
10 RETURN
C ****************************
END
FUNCTION PARDIS (IPRTN, XX, QQ)
C
C Given the parton distribution function in the array U in
C COMMON / PEVLDT / , this routine fetches u(fl, x, q) at any value of
C x and q using Mth-order polinomial (II=0) or rational fraction (II=1)
C interpolation for x. It always uses quadratic polinomial interpolation
C in ln ln (Q/lambda).
C The calling program must ensure that 0 =< x =< 1 ;
C If 0 =< x < Xmin, extrapolation is used and a warning is given
C The calling program must ensure that Alambda < Q ;
C If (Alambda < Q < Qini .or. Qmax < Q),
C extrapolation is used and a warning is given
C
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
Character Msg*80
LOGICAL LSTX
C
PARAMETER (MXX = 1050, MXQ = 25, MXF = 6)
PARAMETER (MXPN = MXF * 2 + 2)
PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
PARAMETER (Smll = 1D-9)
C
COMMON / IOUNIT / NIN, NOUT, NWRT
COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX, DEL,IV
COMMON / XYARAY / ZZ(MXX, MXX), ZV(0:MXX)
COMMON / QARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
COMMON / QARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx
COMMON / PEVLDT / UPD(MXPQX), KF, Nelmt
C
Dimension Fq(5), Df(5)
C M determines the order of the polynomial
C II switches between the polint/ratint interpolation routines
C They are fixed for this version
Save
Data M , II / 2, 0 /
Data Iwrn1, Iwarn2, Iwarn3 / 3*0 /
X = XX
Q = QQ
Md = M / 2
Amd= Md
C
IF (Q .LT. QINI-Smll) THEN
Msg = 'Q less than QINI in PARDIS call; Q SET = QINI.'
CALL WARNR (IWRN1, NWRT, Msg, 'Q', Q, QINI, QMAX, 1)
Q = QINI
ElseIF (Q .GT. QMAX) THEN
Msg = 'Q greater than QMAX in PARDIS call; '
Msg = Msg // 'Extrapolation will be used.'
CALL WARNR(IWRN2, NWRT, Msg, 'Q', Q, QINI, QMAX, 1)
EndIf
C Find lower end of interval containing X
JL = -1
JU = Nx+1
11 If (JU-JL .GT. 1) Then
JM = (JU+JL) / 2
If (X .GT. XV(JM)) Then
JL = JM
Else
JU = JM
Endif
Goto 11
Endif
Jx = JL - (M-1)/2
If (Jx .LT. 0) Then
Jx = 0
Elseif (Jx .GT. Nx-M) Then
Jx = Nx - M
Endif
C Find the interval where Q lies
JL = -1
JU = NT+1
12 If (JU-JL .GT. 1) Then
JM = (JU+JL) / 2
If (Q .GT. QV(JM)) Then
JL = JM
Else
JU = JM
Endif
Goto 12
Endif
Jq = JL - (M-1)/2
If (Jq .LT. 0) Then
Jq = 0
Elseif (Jq .GT. Nt-M) Then
Jq = Nt - M
Endif
SS = LOG (Q/AL)
TT = LOG (SS)
C Find the off-set in the linear array Upd
JFL = IPRTN + NfMx
J0 = (JFL * (NT+1) + Jq) * (NX+1) + Jx + M/2
C
C Nt -| ..........................
C | ..........................
C Jq+M -| .....o ......o............ Iq=M+1
C | .........X................
C Jq -| .(J0)o ......o ........... Iq=1
C | ...... ...................
C 0 --------|-------|-----------|
C 0 Jx Jx+M Nx
Do 21 Iq = 1, M+1
J1 = J0 + (Nx+1)*(Iq-1) + 1
Fq(Iq) = Upd(J1)
21 Continue
C Interpolate in LnLnQ
Call Polint (TV(Jq), Fq(1), M+1, TT, Ftmp, Ddf)
PARDIS = Ftmp
C
RETURN
C ****************************
END
SUBROUTINE EVOLVE (NPO,FINI, IRET)
C ====================================================================
C SUBROUTINE EVOLVE (NPO,FINI, IRET)
C
C Input argument: FINI is a function
C
C FINI (LPARTN, X)
C
C where LPARTN = -6, ... 6 labels the parton flavor:
C t-bar(-6), b-bar(-5), ... gluon(0), u(1), ... t(6) res.
C
C Output parameter:
C
C Iret = 0 : normal execution
C
C NPO : 1 = evolution of H and \tilde H
C 2 = evolution of E and \tilde E
C ----------------------------
C
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
LOGICAL LSTX
C
PARAMETER (MXX = 1050, MXQ = 25, MXF = 6)
PARAMETER (MXPN = MXF * 2 + 2)
PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
C
COMMON / IOUNIT / NIN, NOUT, NWRT
COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX, DEL,IV
COMMON / QARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
COMMON / QARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx
COMMON / PEVLDT / UPD(MXPQX), Kf, Nelmt
C
COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
COMMON / VARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2)
C
DIMENSION QRKP(MXF)
DIMENSION JI(-MXF : MXF+1)
C
EXTERNAL NSRHSP, NSRHSM, FINI
DATA D0 / 0.0 /
C
11 IRET = 0
C Set up number of "valence quarks"
IF (IHDN .LE. 4) THEN
MXVAL = 2
ElseIF (IHDN .LE. 6) THEN
MXVAL = 3
EndIf
C Set up X-mesh points and common parameters
IF (.NOT. LSTX) CALL XARRAY
DLX = 1D0 / NX
J10 = NX / 10
C Set up Q-related arrays
CALL PARPDF (2, 'ALAM', AL, IR)
C Nini and NfMx determined by Qini and Qmax in Qarray
CALL QARRAY (NINI)
C Flavor # NFSN is for the singlet quark combination
NFSN = NFMX + 1
C Kf is the range of the flavor index for the Upd data array
C = quark flavors + 1 for singlet combination and 1 for gluon
KF = 2 * NFMX + 2
C Total number data points to be stored in Upd
Nelmt = Kf * (Nt+1) * (Nx+1)
C ----------------------
C Various initiations:
C
C Offset, or Starting address -1, of Pdf(Iflv) in Upd for the first stage.
C Will be updated each turn
DO 101 IFLV = -NFMX, NFMX+1
JFL = NFMX + IFLV
JI(IFLV) = JFL * (NT+1) * (NX+1)
101 CONTINUE
C Define initial distributions for the evolution.
C Input functions are defined on (1:Nx).
C
C Input gluon and quark distributions are inserted
C in the Upd array directly; Output of each stage of
C evolution serve as the input for the next stage.
C
3 DO 31 IZ = 1, NX
C Gluon: (Position Ji(0)+1 is reserved for x = 0)
C
UPD(JI(0)+IZ+1) = FINI (NPO,0, XV(IZ))
C Singlet quark
C Input singlet quark distribution must be initiated for
C each stage separately.
C
UPD(JI(NFSN)+IZ+1) = 0
C If no quark, bypass filling of quark arrays
IF (NFMX .EQ. 0) GOTO 31
C
DO 331 IFLV = 1, NINI
A = FINI (NPO, IFLV, XV(IZ))
B = FINI (NPO,-IFLV, XV(IZ))
IF (XV(IZ).GT.DEL) THEN
QRKP (IFLV) = A + B
ELSE
QRKP(IFLV) = A
ENDIF
C Acculumate singlet distribution
UPD(JI(NFSN)+IZ+1) = UPD(JI(NFSN)+IZ+1) + QRKP (IFLV)
C Initialize the "minus" non-singlet com-
C bination (Q - Q-bar) in area for Q-bar
IF (XV(IZ).GT.DEL) THEN
UPD(JI(-IFLV)+IZ+1) = A - B
ELSE
UPD(JI(-IFLV)+IZ+1) = B
ENDIF
331 CONTINUE
C The "plus" non-singlet combination is initialized
C in the array area for the Quark of the same flavor
DO 332 IFLV = 1, NINI
UPD(JI(IFLV)+IZ+1) = QRKP(IFLV) - UPD(JI(NFSN)+IZ+1)/NINI
332 CONTINUE
C
31 CONTINUE
C -----------------------
C Start of the Q2- Evolution:
C
C Outer loop is by the effective number of flvr
DO 21 NEFF = NINI, NFMX
C Set up 1st and 2nd order kernel functions
IF (IKNL .EQ. 2 .OR. IKNL.EQ.-2) CALL STUPKL (NEFF)
C Singlet Calculation
ICNT = NEFF - NINI + 1
C Skip if new quark mass = old
IF (NTN(ICNT) .EQ. 0) GOTO 21
C Otherwise, recall iteration parameters
NITR = NTN (ICNT)
DT = DTN (ICNT)
TIN = TLN (ICNT)
C Perform evolution
CALL SNEVL (NPO,XV,IV,DEL,IKNL, NX, NITR, JT, DT, TIN, NEFF,
> UPD(JI(NFSN)+2), UPD(JI(0)+2), UPD(JI(NFSN)+1), UPD(JI(0)+1))
C
C Non-singlet sector, one flavor a time.
C
C Skip this section if quark flavor = 0.
IF (NEFF .EQ. 0) GOTO 88
5 DO 333 IFLV = 1, NEFF
C First evolve the (Q+Q-bar) "PLUS NS" part
CALL NSEVL (NPO,0,NSRHSP,XV,IV,DEL,IKNL, NX, NITR, JT, DT,
> TIN, NEFF, UPD(JI( IFLV)+2), UPD(JI( IFLV)+1))
C
C The (Q-Q-bar) "MINUS NS" evolution is needed
C only for those flavors with valence
IF (IFLV .LE. MXVAL)
> CALL NSEVL (NPO,1,NSRHSM,XV,IV,DEL,IKNL, NX, NITR, JT, DT,
> TIN, NEFF, UPD(JI(-IFLV)+2), UPD(JI(-IFLV)+1))
C
C To obtain the real quark distribution functions,
C combine the singlet piece to the two non-singlet pieces.
C Enforce positivity conditions also at this stage.
C
DO 55 IS = 0, NITR
DO 56 IX = 0, NX
TP = UPD (IS*(NX+1) + IX + 1 + JI( IFLV))
TS = UPD (IS*(NX+1) + IX + 1 + JI( NFSN)) / NEFF
TP = TP + TS
IF (IKNL .GT. 0 .AND. XV(IX).GT.DEL) TP = TP
IF (IFLV .LE. MXVAL) THEN
TM = UPD (IS*(NX+1) + IX + 1 + JI(-IFLV))
IF (IKNL .GT. 0 .AND. XV(IX).GT.DEL) THEN
TM = TM
c TP = MAX (TP, TM)
TP = TP
EndIf
Else
TM = 0.
EndIf
IF (XV(IX).GT.DEL) THEN
UPD (JI( IFLV) + IS*(NX+1) + IX + 1) = (TP + TM)/2.
UPD (JI(-IFLV) + IS*(NX+1) + IX + 1) = (TP - TM)/2.
ELSE
UPD (JI( IFLV) + IS*(NX+1) + IX + 1) = TP
UPD (JI(-IFLV) + IS*(NX+1) + IX + 1) = TM
ENDIF
56 CONTINUE
IF (IS.EQ.0) THEN
GOTO 55
ELSE
ENDIF
55 CONTINUE
333 CONTINUE
C Heavy quarks above current threshold have zero distribution
C
DO 334 IFLV = NEFF + 1, NFMX
DO 57 IS = 0, NITR
DO 58 IX = 0, NX
UPD(JI( IFLV) + IS*(NX+1) + IX + 1) = 0
UPD(JI(-IFLV) + IS*(NX+1) + IX + 1) = 0
58 CONTINUE
57 CONTINUE
334 CONTINUE
88 CONTINUE
C ------------------------------
C Define initial parameters for next stage
IF (NFMX .EQ. NEFF) GOTO 21
C New Offsets
DO 335 IFLV = -NFMX, NFMX+1
JI(IFLV) = JI(IFLV) + NITR * (NX+1)
335 CONTINUE
C New distributions:
C gluon input functions are in place;
C only non-singlet and singlet quark needs re-initiation
C
C Calculate initial heavy quark distribution due to
C change of renormalization scheme across threshold:
CALL HQRK (NX, TT, NEFF+1, UPD(JI(0)+2), UPD(JI(NEFF+1)+2))
C
DO 32 IZ = 1, NX
QRKP (NEFF+1) = 2. * UPD(JI( NEFF+1) + IZ + 1)
C New Singlet piece
UPD (JI(NFSN)+IZ+1) = UPD (JI(NFSN)+IZ+1) + QRKP (NEFF+1)
VS00 = UPD (JI(NFSN)+IZ+1) / (NEFF+1)
C
C "plus" non-singlet for the new flavor
UPD(JI( NEFF+1) + IZ + 1) = QRKP(NEFF+1) - VS00
C
C Calculate the non-singulet parts of the other quark distr.
C
C Change from the output of last stage of calcu-
C lation due to two sources of change in Vs/Neff:
C change of Neff and addition of new quark distr.
DO 321 IFL = 1, NEFF
A = UPD(JI( IFL)+IZ+1)
B = UPD(JI(-IFL)+IZ+1)
IF (XV(IZ).GT.DEL) THEN
QRKP(IFL) = A + B
ELSE
QRKP(IFL) = A
ENDIF
C "plus" non-singlet for flavor IFL
UPD(JI( IFL)+IZ+1) = QRKP(IFL) - VS00
C "minus" non-singlet for flavors with valence
IF (XV(IZ).GT.DEL) THEN
IF (IFL .LE. MXVAL) UPD(JI(-IFL)+IZ+1) = A - B
ELSE
IF (IFL .LE. MXVAL) UPD(JI(-IFL)+IZ+1) = B
ENDIF
321 CONTINUE
C
32 CONTINUE
C Return of Q-2 evolution loop to the next stage
21 CONTINUE
C Conclusion of the full calculation
Return
C **********************
End
SUBROUTINE NSEVL (NPO,NSET,RHS,XT,IV,DEL,
> IKNL,NX,NT,JT,DT,TIN,NEFF,U0,UN)
C
C IKNL determines to which order in Alpha is the
C calculation to be carried out.
C
C Given the non-singlet parton distribution function U0 at some initial
C QIN (Tt= 0) in the x interval (0, 1) covered by the array Iz = 1, NX,
C this routine calculates the evoluted function U at Nt values of Tt at
C intervals of Dt by numerically integrating the A-P equation using the
C non-singlet kernel.
C
C Un(Ix, Tt) = Y(x,es) at the sites Ix=0,..,Nx (x = Ix * Dz);
C Un(0, Tt) is obtained by quadratic extrapolation from Ix = 1, 2, 3
C for each Tt rather then by evolution because of possible
C singular behavior at x = 0.
C
C Data is stored at Tt = Is * Dt, Is = 0, ... , Nt.
C The function at Is = 0 is the input distribution.
C
C Numerical iteration is performed with finer grain Ddt = Dt/Jt.
C
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
PARAMETER (MXX = 1050, MXQ = 25, MXF = 6)
PARAMETER (MXPN = MXF * 2 + 2)
PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
C
COMMON / IOUNIT / NIN, NOUT, NWRT
Common / PdCntrl/ LPrt, LDbg
COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
C
DIMENSION U0(NX), UN(0:NX, 0:NT)
DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX)
DIMENSION XT(0:MXX)
external rhs
C
If (Ldbg .Eq. 1) Then
Write (Nwrt, '(A)') ' Non-Singlet:'
N5 = Nx / 5 + 1
Endif
DDT = DT / JT
C
IF (NX .GT. MXX) THEN
WRITE (NOUT,*) 'Nx =', NX, ' greater than Max # of pts in NSEVL.'
STOP 'Program stopped in NSEVL'
EndIf
C
C Compute an effective first order lamda
C to be used in checking of moment evl.
TMD = TIN + DT * NT / 2.
AMU = EXP(TMD)
TEM = 6./ (33.- 2.* NEFF) / ALPI(AMU)
TLAM = TMD - TEM
C
C Fill first rows of output array, for initial value of Q
C
DO 9 IX = 1, NX
UN(IX, 0) = U0(IX)
9 CONTINUE
UN(0, 0) = 3D0*U0(1) - 3D0*U0(2) - U0(1)
C
C Initiation
C
TT = TIN
DO 10 IZ = 1, NX
Y0(IZ) = U0(IZ)
10 CONTINUE
C
C loop in the Tt variable
C
DO 20 IS = 1, NT
C
C fine- grained iteration
C
IPO = JT
DO 202 JS = 1, JT
C
C Irnd is the counter of the Q-iteration
C
IRND = (IS-1) * JT + JS
C
C Use Runge-Katta the first round
C
IF (IRND .EQ. 1) THEN
C
CALL RHS (NPO,TT, Neff, Y0, F0)
C smoothing
CALL SMOOTH (NSET,IKNL,MXX,DEL,IV,XT,F0)
DO 250 IZ = 1, NX
Y0(IZ) = Y0(IZ) + DDT * F0(IZ)
250 CONTINUE
C smoothing
CALL SMOOTH (NSET,IKNL,MXX,DEL,IV,XT,Y0)
C
TT = TT + DDT
C
CALL RHS (NPO,TT, NEFF, Y0, F1)
C smoothing
CALL SMOOTH (NSET,IKNL,MXX,DEL,IV,XT,F1)
DO 251 IZ = 1, NX
Y1(IZ) = U0(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0
251 CONTINUE
c smoothing
CALL SMOOTH (NSET,IKNL,MXX,DEL,IV,XT,Y1)
C
C What follows is a combination of the 2-step method
C plus the Adams Predictor-Corrector Algorithm
C
Else
C
CALL RHS (NPO,TT, NEFF, Y1, F1)
c smoothing
CALL SMOOTH (NSET,IKNL,MXX,DEL,IV,XT,F1)
C
C Predictor
C
DO 252 IZ = 1, NX
YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0
252 CONTINUE
c smoothing
CALL SMOOTH (NSET,IKNL,MXX,DEL,IV,XT,YP)
C
C Increment of Tt at this place is part of the formalism
C
TT = TT + DDT
C
CALL RHS (NPO,TT, NEFF, YP, FP)
C smoothing
CALL SMOOTH (NSET,IKNL,MXX,DEL,IV,XT,FP)
C
C Corrector
C
DO 253 IZ = 1, NX
Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0
F0(IZ) = F1(IZ)
253 CONTINUE
c smoothing
CALL SMOOTH (NSET,IKNL,MXX,DEL,IV,XT,Y1)
EndIf
C
202 CONTINUE
C
C Fill output array
C
DO 260 IZ = 1, NX
UN (IZ, IS) = Y1(IZ)
260 CONTINUE
C
C The value of the function at x=0 is obtained by extrapolation
C
UN(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3)
C
C Print out for Debugging
C
If (LDbg .Eq. 1) Then
Write (Nwrt, '(A, 5(1pE12.3))') ' :', (Un(Iz,Is), Iz=1,Nx,N5)
Endif
C
JT = IPO
20 CONTINUE
C
RETURN
C ****************************
END
SUBROUTINE NSRHSM (NPO,TT, NEFF, FI, FO)
C
C Subroutine to compute the Right-Side of the Altarelli-Parisi Equation
C This copy applies to the "NS-minus" piece -- (Qrk - Qrk-bar)
C
C See comments in NSRHSP for details
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
LOGICAL LSTX
C
PARAMETER (MXX = 1050, MXQ = 25, MXF = 6)
PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
PARAMETER (MXQX = MXX*MXQ)
PARAMETER (PI = 3.141592653589793, PI2 = PI ** 2,
>Z3 = 2.404113806319188)
PARAMETER (D0 = 0.0, D1 = 1.0)
COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX, DEL, IV
COMMON / XYARAY / ZZ(MXX, MXX), ZV(0:MXX)
COMMON / KRNL01 / AGG2 (0:MXX), ANSP (0:MXX), ANSM (0:MXX)
COMMON / KRN1ST / FF1(0:MXX,0:MXX), GG1(0:MXX,0:MXX),
> FG2(0:MXX,0:MXX),GF2(0:MXX,0:MXX),GG2(0:MXX,0:MXX),
> PNSP(0:MXX,0:MXX), PNSM(0:MXX,0:MXX), SFF2(0:MXX,0:MXX),
> SFF2BLK(0:MXX,0:MXX), FG2BLK(0:MXX,0:MXX), GF2BLK(0:MXX,0:MXX),
> GG2BLK(0:MXX,0:MXX), PNSPBLK(0:MXX,0:MXX), PNSMBLK(0:MXX,0:MXX),
> GG2BLC1(0:MXX,0:MXX), PNSPBLC1(0:MXX,0:MXX),FF1BLX(0:MXX,0:MXX),
> GG2BLCX(0:MXX,0:MXX), PNSPBLCX(0:MXX,0:MXX),GG1BLX(0:MXX,0:MXX),
> PNSMBLC1(0:MXX,0:MXX), PNSMBLCX(0:MXX,0:MXX),FF1BL1(0:MXX,0:MXX),
> GG1BL1(0:MXX,0:MXX),SFF2BLKD(0:MXX,0:MXX), FG2BLKD(0:MXX,0:MXX),
> GF2BLKD(0:MXX,0:MXX), GG2BLKD(0:MXX,0:MXX),PNSPBLKD(0:MXX,0:MXX),
> PNSMBLKD(0:MXX,0:MXX),FG2BLK1(0:MXX,0:MXX),GF2BLK1(0:MXX,0:MXX),
> SFF2BLK1(0:MXX,0:MXX),PNSPC(0:MXX,0:MXX),PNSMC(0:MXX,0:MXX),
> GG2C(0:MXX,0:MXX), PNSPDGBL(0:MXX,0:MXX), PNSMDGBL(0:MXX,0:MXX),
> GG2DGBL(0:MXX,0:MXX)
COMMON / REGFUNC / GB1(MXX), GBX(MXX), GB21(MXX), GB510(MXX),
> GB2X(MXX), GB31(MXX), GB41(MXX), GB5X(MXX), GB51(MXX),
> GB3X(MXX), GB4X(MXX), GB210(MXX), GB10(MXX), GB310(MXX),
> G1(MXX), G2(MXX), G3(MXX), G4(MXX), G5(MXX), GB410(MXX)
COMMON / DIFFD / DIFFDEL
c COMMON / SECDERIV / DGBX(3,MXX), DGB2X(3,MXX), DGB1(3,MXX),
c > DGB21(3,MXX),
c > DGB01(3,MXX), DGB41(3,MXX), DGB10(3,MXX), DGB210(3,MXX),
c > DG1(3,MXX), DG2(3,MXX)
COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx
DIMENSION FI(NX), FO(NX), FN(MXX), FFN(3*MXX),FQ(MXX)
DIMENSION W0(MXX), W1(MXX), WH(MXX), W2(MXX), W01(MXX)
DIMENSION W3(MXX), W5(MXX), W22(MXX), FFI(3*MXX), WH1(MXX)
DIMENSION X2(10), X1(10), FX(10), FX1(10), AF(MXX), AF1(MXX),
> AF2(MXX), TQQ1(MXX), TQQX(MXX),TQQ11(MXX),TQQD(MXX)
DIMENSION TEMPI1(0:MXX),TEMPI2(0:MXX)
EXTERNAL SMPSNA,SMPNOL,SMPNOR
SAVE
DATA AERR, RERR / 0.00001, 0.0002/
C
S = EXP(TT)
Q = AL * EXP (S)
CPL = ALPI(Q)
CPL2= CPL ** 2 / 2. * S
CPL = CPL * S
C
C Set up Quark pdf for integration in BL region!
C
DZ = 1./(NX-1)
FN(NX) = 0.0
FQ(NX) = 0.0
DO 229 IY = 1, NX-1
FQ(IY) = FI(IY)*XA(IY,1)
TMX = ABS(XV(IY)-DEL/2.)
IF (TMX.GT.1E-10) then
FN(IY) = FI(IY)
ELSE
IF (IKNL.LT.0) then
FN(IY) = 0.0
ELSE
FN(IY) = FI(IY)
ENDIF
ENDIF
229 CONTINUE
C
C LO convolution of kernel with pdf! First unpol. then pol. case!
C
CALL NEWARRAY (NX, DEL, FQ, FFI)
CALL NEWARRAY (NX, DEL, FN, FFN)
CALL INTEGR (NX, 0, FQ, W0, IR1)
CALL INTEGR (NX, 0, FN, W01, IR12)
CALL NINTEGR (NX, DEL, 1, FFI, FQ, W1, IR11)
CALL NINTEGR (NX, DEL, 15, FFN, FN, W2, IR01)
CALL NINTEGR (NX, DEL, 8, FFI, FQ, W22, IR02)
CALL NINTEGR (NX, DEL, 7, FFI, FQ, W3, IR22)
CALL NINTEGR (NX, DEL, 12, FFN, FN, W5, IR3)
CALL NINTEGR (NX, DEL, 16, FFN, FN, AF, IR4)
CALL NINTEGR (NX, DEL, 20, FFN, FN, AF2, IR4)
CALL HINTEGN (2,NX,IV,DEL,FFI,FQ, WH)
CALL HINTEGN (1,NX,IV,DEL,FFN,FN, WH1)
C
DIM1 = 1D0/DEL
DO 230 IZ = 2, NX
C
C First treatment of ERBL region
C
IF (XA(IZ,1) .LE. DEL) THEN
AA1 = 1D0 - XA(IZ,1)*DIM1
C
C Special treatment of the point y=del and y= 0
C
IF (IZ.EQ.IV) THEN
C
C For H and \tilde H use normal extrapolation and then symmetry properties
C for E and \tilde E set those points to zero!