source: LMDZ4/trunk/libf/phylmd/cv3_inip.F @ 1085

Last change on this file since 1085 was 987, checked in by Laurent Fairhead, 16 years ago

Du nettoyage sur le parallelisme, inclusion de nouvelles interfaces pour OPA9
et ORCHIDEE YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.7 KB
RevLine 
[963]1        SUBROUTINE cv3_inip()
[879]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"
12c
[963]13c      INTEGER iflag_mix
[879]14c
15c --   Mixing probability distribution functions
16c
17      real Qcoef1,Qcoef2,QFF,QFFF,Qmix,Rmix,Qmix1,Rmix1,Qmix2,Rmix2,F
18      Qcoef1(F) = tanh(F/gammas)
19      Qcoef2(F) = ( tanh(F/gammas) + gammas *
20     $            log(cosh((1.- F)/gammas)/cosh(F/gammas)))
21      QFF(F) = Max(Min(F,1.),0.)
22      QFFF(F) = Min(QFF(F),scut)
23      Qmix1(F) = ( tanh((QFF(F) - Fmax)/gammas)+Qcoef1max )/
24     $           Qcoef2max
25      Rmix1(F) = ( gammas*log(cosh((QFF(F)-Fmax)/gammas))
26     1             +QFF(F)*Qcoef1max ) / Qcoef2max
27      Qmix2(F) = -Log(1.-QFFF(F))/scut
28      Rmix2(F) = (QFFF(F)+(1.-QFF(F))*Log(1.-QFFF(F)))/scut
29      Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
30      Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
31C
32ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
33c
34C
35C===========================================================================
36C       READ IN PARAMETERS FOR THE MIXING DISTRIBUTION
37C       AND PASS THESE THROUGH A COMMON BLOCK TO SUBROUTINE CONVECT etc.
38C       (Written by V.T.J. Phillips, 20-30/Jan/99)
39C===========================================================================
40C
41C   line 1:  a flag (0 or 1) to decide whether P(F) = 1 or the general P(F) is to be
42C         used, followed by SCUT, which is the cut-off value of F in CONVECT
43C   line 2:  blank
44C   line 3:  the coefficients for the linear combination of P(F)s to
45C                 make the general P(F)
46C   line 4:  blank
47C   line 5:  gammas, Fmax for the cosh^2 component of P(F)
48C   line 6:  blank
49C   line 7:  alphas for the 1st irrational P(F)
50C   line 8:  blank
51C   line 9:  betas  for the 2nd irrational P(F)
52C
53
[987]54cc$$$        open(57,file='parameter_mix.data')
55cc$$$
56cc$$$        read(57,*) iflag_mix, scut
57cc$$$        read(57,*)
58cc$$$        if(iflag_mix .gt. 0) then
59cc$$$         read(57,*) qqa1, qqa2
60cc$$$              read(57,*)
61cc$$$              read(57,*) gammas, Fmax
62cc$$$              read(57,*)
63cc$$$              read(57,*) alphas
64cc$$$         endif
65cc$$$    close(57)
[879]66
67c
68      if(iflag_mix .gt. 0) then
69c
70c--      Normalize Pdf weights
71c
72        sumcoef=qqa1+qqa2
73        qqa1=qqa1/sumcoef
74        qqa2=qqa2/sumcoef
75c
76        Qcoef1max = Qcoef1(Fmax)
77        Qcoef2max = Qcoef2(Fmax)
78c
79        sigma = 0.
80        aire=0.0
81        pdf=0.0
82        mu=0.0
83        df = 0.0001
84c
85c        do ff = 0.0 + df, 1.0 - 2.*df, df
86         ff=df
87         dowhile ( ff .le. 1.0 - 2.*df )
88              pdf = (Qmix(ff+df) -  Qmix(ff)) * (1.-ff) / df
89              aire=aire+(Qmix(ff+df) - Qmix(ff)) * (1.-ff)
90              mu = mu + pdf * ff * df
91              write(*,*) pdf,  Qmix(ff), aire, ff
92         ff=ff+df
93         enddo
94c
95c         do ff=0.0+df,1.0 - 2.*df,df
96          ff=df
97          dowhile ( ff .le. 1.0 - 2.*df )
98              pdf = (Qmix(ff+df)- Qmix(ff)) * (1.-ff) / df
99              sigma = sigma+pdf*(ff - mu)*(ff - mu)*df
100         ff=ff+df
101         enddo
102         sigma = sqrt(sigma)
103c
104        if (abs(aire-1.0) .gt. 0.02) then
105            print *,'WARNING:: AREA OF MIXING PDF IS::', aire
106            stop
107        else
108            print *,'Area, mean & std deviation are ::', aire,mu,sigma
109        endif
110      endif     !  (iflag_mix .gt. 0)
111
112      RETURN
113      END
Note: See TracBrowser for help on using the repository browser.