9 common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
13 common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
21 IF(mvar.EQ.3.OR.mvar.EQ.4)
THEN
24 ELSEIF(mvar.EQ.5.OR.mvar.EQ.6)
THEN
28 IF(
mint(43).EQ.1.AND.(iset(isub).EQ.1.OR.iset(isub).EQ.2))
THEN
30 ELSEIF(mvar.EQ.1)
THEN
31 tau=taumin*(taumax/taumin)**vvar
32 ELSEIF(mvar.EQ.2)
THEN
33 tau=taumax*taumin/(taumin+(taumax-taumin)*vvar)
34 ELSEIF(mvar.EQ.3.OR.mvar.EQ.5)
THEN
35 ratgen=(taure+taumax)/(taure+taumin)*taumin/taumax
36 tau=taure*taumin/((taure+taumin)*ratgen**vvar-taumin)
38 aupp=atan((taumax-taure)/gamre)
39 alow=atan((taumin-taure)/gamre)
40 tau=taure+gamre*tan(alow+(aupp-alow)*vvar)
45 ELSEIF(ivar.EQ.2)
THEN
48 IF(
mint(43).EQ.1)
THEN
50 ELSEIF(
mint(43).EQ.2)
THEN
51 IF(iset(isub).LE.2) yst=-0.5*
log(
vint(21))
52 IF(iset(isub).GE.3) yst=-0.5*
log(
vint(26))
53 ELSEIF(
mint(43).EQ.3)
THEN
54 IF(iset(isub).LE.2) yst=0.5*
log(
vint(21))
55 IF(iset(isub).GE.3) yst=0.5*
log(
vint(26))
56 ELSEIF(mvar.EQ.1)
THEN
57 yst=ystmin+(ystmax-ystmin)*sqrt(vvar)
58 ELSEIF(mvar.EQ.2)
THEN
59 yst=ystmax-(ystmax-ystmin)*sqrt(1.-vvar)
61 aupp=atan(exp(ystmax))
62 alow=atan(exp(ystmin))
63 yst=
log(tan(alow+(aupp-alow)*vvar))
65 vint(22)=min(ystmax,
max(ystmin,yst))
68 ELSEIF(ivar.EQ.3)
THEN
80 IF(aneg.GT.0..AND.vvar*(aneg+apos).LE.aneg)
THEN
81 vctn=vvar*(aneg+apos)/aneg
82 cth=ctnmin+(ctnmax-ctnmin)*vctn
84 vctp=(vvar*(aneg+apos)-aneg)/apos
85 cth=ctpmin+(ctpmax-ctpmin)*vctp
87 ELSEIF(mvar.EQ.2)
THEN
88 rmnmin=
max(rm34,rsqm-ctnmin)
89 rmnmax=
max(rm34,rsqm-ctnmax)
90 rmpmin=
max(rm34,rsqm-ctpmin)
91 rmpmax=
max(rm34,rsqm-ctpmax)
92 aneg=
log(rmnmin/rmnmax)
93 apos=
log(rmpmin/rmpmax)
94 IF(aneg.GT.0..AND.vvar*(aneg+apos).LE.aneg)
THEN
95 vctn=vvar*(aneg+apos)/aneg
96 cth=rsqm-rmnmin*(rmnmax/rmnmin)**vctn
98 vctp=(vvar*(aneg+apos)-aneg)/apos
99 cth=rsqm-rmpmin*(rmpmax/rmpmin)**vctp
101 ELSEIF(mvar.EQ.3)
THEN
102 rmnmin=
max(rm34,rsqm+ctnmin)
103 rmnmax=
max(rm34,rsqm+ctnmax)
104 rmpmin=
max(rm34,rsqm+ctpmin)
105 rmpmax=
max(rm34,rsqm+ctpmax)
106 aneg=
log(rmnmax/rmnmin)
107 apos=
log(rmpmax/rmpmin)
108 IF(aneg.GT.0..AND.vvar*(aneg+apos).LE.aneg)
THEN
109 vctn=vvar*(aneg+apos)/aneg
110 cth=rmnmin*(rmnmax/rmnmin)**vctn-rsqm
112 vctp=(vvar*(aneg+apos)-aneg)/apos
113 cth=rmpmin*(rmpmax/rmpmin)**vctp-rsqm
115 ELSEIF(mvar.EQ.4)
THEN
116 rmnmin=
max(rm34,rsqm-ctnmin)
117 rmnmax=
max(rm34,rsqm-ctnmax)
118 rmpmin=
max(rm34,rsqm-ctpmin)
119 rmpmax=
max(rm34,rsqm-ctpmax)
120 aneg=1./rmnmax-1./rmnmin
121 apos=1./rmpmax-1./rmpmin
122 IF(aneg.GT.0..AND.vvar*(aneg+apos).LE.aneg)
THEN
123 vctn=vvar*(aneg+apos)/aneg
124 cth=rsqm-1./(1./rmnmin+aneg*vctn)
126 vctp=(vvar*(aneg+apos)-aneg)/apos
127 cth=rsqm-1./(1./rmpmin+apos*vctp)
129 ELSEIF(mvar.EQ.5)
THEN
130 rmnmin=
max(rm34,rsqm+ctnmin)
131 rmnmax=
max(rm34,rsqm+ctnmax)
132 rmpmin=
max(rm34,rsqm+ctpmin)
133 rmpmax=
max(rm34,rsqm+ctpmax)
134 aneg=1./rmnmin-1./rmnmax
135 apos=1./rmpmin-1./rmpmax
136 IF(aneg.GT.0..AND.vvar*(aneg+apos).LE.aneg)
THEN
137 vctn=vvar*(aneg+apos)/aneg
138 cth=1./(1./rmnmin-aneg*vctn)-rsqm
140 vctp=(vvar*(aneg+apos)-aneg)/apos
141 cth=1./(1./rmpmin-apos*vctp)-rsqm
144 IF(cth.LT.0.) cth=min(ctnmax,
max(ctnmin,cth))
145 IF(cth.GT.0.) cth=min(ctpmax,
max(ctpmin,cth))
149 ELSEIF(ivar.EQ.4)
THEN
153 IF(
mint(43).EQ.1)
THEN
155 ELSEIF(mvar.EQ.1)
THEN
156 taup=taupmn*(taupmx/taupmn)**vvar
158 aupp=(1.-
tau/taupmx)**4
159 alow=(1.-
tau/taupmn)**4
160 taup=
tau/(1.-(alow+(aupp-alow)*vvar)**0.25)
162 vint(26)=min(taupmx,
max(taupmn,taup))