1 |
|
|
SUBROUTINE cv3_inicp() |
2 |
|
|
|
3 |
|
|
! ************************************************************** |
4 |
|
|
! * |
5 |
|
|
! CV3_INIP Lecture des choix de lois de probabilit� de m�lange* |
6 |
|
|
! et calcul de leurs coefficients normalis�s. * |
7 |
|
|
! * |
8 |
|
|
! written by : Jean-Yves Grandpeix, 06/06/2006, 19.39.27 * |
9 |
|
|
! modified by : * |
10 |
|
|
! ************************************************************** |
11 |
|
|
IMPLICIT NONE |
12 |
|
|
include "YOMCST2.h" |
13 |
|
|
|
14 |
|
|
INTEGER iflag_clos |
15 |
|
|
CHARACTER (LEN=20) :: modname = 'cv3_inicp' |
16 |
|
|
CHARACTER (LEN=80) :: abort_message |
17 |
|
|
|
18 |
|
|
! -- Mixing probability distribution functions |
19 |
|
|
|
20 |
|
|
REAL qcoef1, qcoef2, qff, qfff, qmix, rmix, qmix1, rmix1, qmix2, rmix2, f |
21 |
|
|
REAL :: sumcoef,sigma,aire,pdf,mu,df,ff |
22 |
|
|
|
23 |
|
|
qcoef1(f) = tanh(f/gammas) |
24 |
|
|
qcoef2(f) = (tanh(f/gammas)+gammas*log(cosh((1.-f)/gammas)/cosh(f/gammas))) |
25 |
|
|
qff(f) = max(min(f,1.), 0.) |
26 |
|
|
qfff(f) = min(qff(f), scut) |
27 |
|
|
qmix1(f) = (tanh((qff(f)-fmax)/gammas)+qcoef1max)/qcoef2max |
28 |
|
|
rmix1(f) = (gammas*log(cosh((qff(f)-fmax)/gammas))+qff(f)*qcoef1max)/ & |
29 |
|
|
qcoef2max |
30 |
|
|
qmix2(f) = -log(1.-qfff(f))/scut |
31 |
|
|
rmix2(f) = (qfff(f)+(1.-qff(f))*log(1.-qfff(f)))/scut |
32 |
|
|
qmix(f) = qqa1*qmix1(f) + qqa2*qmix2(f) |
33 |
|
|
rmix(f) = qqa1*rmix1(f) + qqa2*rmix2(f) |
34 |
|
|
|
35 |
|
|
! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
36 |
|
|
|
37 |
|
|
|
38 |
|
|
! =========================================================================== |
39 |
|
|
! READ IN PARAMETERS FOR THE MIXING DISTRIBUTION |
40 |
|
|
! AND PASS THESE THROUGH A COMMON BLOCK TO SUBROUTINE CONVECT etc. |
41 |
|
|
! (Written by V.T.J. Phillips, 20-30/Jan/99) |
42 |
|
|
! =========================================================================== |
43 |
|
|
|
44 |
|
|
! line 1: a flag (0 or 1) to decide whether P(F) = 1 or the general P(F) |
45 |
|
|
! is to be |
46 |
|
|
! used, followed by SCUT, which is the cut-off value of F in CONVECT |
47 |
|
|
! line 2: blank |
48 |
|
|
! line 3: the coefficients for the linear combination of P(F)s to |
49 |
|
|
! make the general P(F) |
50 |
|
|
! line 4: blank |
51 |
|
|
! line 5: gammas, Fmax for the cosh^2 component of P(F) |
52 |
|
|
! line 6: blank |
53 |
|
|
! line 7: alphas for the 1st irrational P(F) |
54 |
|
|
! line 8: blank |
55 |
|
|
! line 9: betas for the 2nd irrational P(F) |
56 |
|
|
|
57 |
|
|
|
58 |
|
|
! open(57,file='parameter_mix.data') |
59 |
|
|
|
60 |
|
|
! read(57,*) iflag_clos |
61 |
|
|
! read(57,*) iflag_mix, scut |
62 |
|
|
! read(57,*) |
63 |
|
|
! if(iflag_mix .gt. 0) then |
64 |
|
|
! read(57,*) qqa1, qqa2 |
65 |
|
|
! read(57,*) |
66 |
|
|
! read(57,*) gammas, Fmax |
67 |
|
|
! read(57,*) |
68 |
|
|
! read(57,*) alphas |
69 |
|
|
! endif |
70 |
|
|
! close(57) |
71 |
|
|
|
72 |
|
|
IF (iflag_mix>0) THEN |
73 |
|
|
|
74 |
|
|
! -- Normalize Pdf weights |
75 |
|
|
|
76 |
|
|
sumcoef = qqa1 + qqa2 |
77 |
|
|
qqa1 = qqa1/sumcoef |
78 |
|
|
qqa2 = qqa2/sumcoef |
79 |
|
|
|
80 |
|
|
qcoef1max = qcoef1(fmax) |
81 |
|
|
qcoef2max = qcoef2(fmax) |
82 |
|
|
|
83 |
|
|
sigma = 0. |
84 |
|
|
aire = 0.0 |
85 |
|
|
pdf = 0.0 |
86 |
|
|
mu = 0.0 |
87 |
|
|
df = 0.0001 |
88 |
|
|
|
89 |
|
|
! do ff = 0.0 + df, 1.0 - 2.*df, df |
90 |
|
|
ff = df |
91 |
|
|
DO WHILE (ff<=1.0-2.*df) |
92 |
|
|
pdf = (qmix(ff+df)-qmix(ff))*(1.-ff)/df |
93 |
|
|
aire = aire + (qmix(ff+df)-qmix(ff))*(1.-ff) |
94 |
|
|
mu = mu + pdf*ff*df |
95 |
|
|
! c write(*,*) pdf, Qmix(ff), aire, ff |
96 |
|
|
ff = ff + df |
97 |
|
|
END DO |
98 |
|
|
|
99 |
|
|
! do ff=0.0+df,1.0 - 2.*df,df |
100 |
|
|
ff = df |
101 |
|
|
DO WHILE (ff<=1.0-2.*df) |
102 |
|
|
pdf = (qmix(ff+df)-qmix(ff))*(1.-ff)/df |
103 |
|
|
sigma = sigma + pdf*(ff-mu)*(ff-mu)*df |
104 |
|
|
ff = ff + df |
105 |
|
|
END DO |
106 |
|
|
sigma = sqrt(sigma) |
107 |
|
|
|
108 |
|
|
IF (abs(aire-1.0)>0.02) THEN |
109 |
|
|
PRINT *, 'WARNING:: AREA OF MIXING PDF IS::', aire |
110 |
|
|
abort_message = '' |
111 |
|
|
CALL abort_physic(modname, abort_message, 1) |
112 |
|
|
ELSE |
113 |
|
|
PRINT *, 'Area, mean & std deviation are ::', aire, mu, sigma |
114 |
|
|
END IF |
115 |
|
|
END IF ! (iflag_mix .gt. 0) |
116 |
|
|
|
117 |
|
|
RETURN |
118 |
|
|
END SUBROUTINE cv3_inicp |