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