source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/cv3_inip.F90 @ 2283

Last change on this file since 2283 was 1999, checked in by Laurent Fairhead, 11 years ago

Merged trunk changes r1920:1997 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.5 KB
Line 
1SUBROUTINE 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
13  ! INTEGER iflag_mix
14  include 'iniprint.h'
15
16  CHARACTER (LEN=20) :: modname = 'cv3_inip'
17  CHARACTER (LEN=80) :: abort_message
18
19
20  ! --   Mixing probability distribution functions
21
22  REAL qcoef1, qcoef2, qff, qfff, qmix, rmix, qmix1, rmix1, qmix2, rmix2, f
23
24  qcoef1(f) = tanh(f/gammas)
25  qcoef2(f) = (tanh(f/gammas)+gammas*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)/qcoef2max
29  rmix1(f) = (gammas*log(cosh((qff(f)-fmax)/gammas))+qff(f)*qcoef1max)/ &
30    qcoef2max
31  qmix2(f) = -log(1.-qfff(f))/scut
32  rmix2(f) = (qfff(f)+(1.-qff(f))*log(1.-qfff(f)))/scut
33  qmix(f) = qqa1*qmix1(f) + qqa2*qmix2(f)
34  rmix(f) = qqa1*rmix1(f) + qqa2*rmix2(f)
35
36  ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
37
38
39  ! ===========================================================================
40  ! READ IN PARAMETERS FOR THE MIXING DISTRIBUTION
41  ! AND PASS THESE THROUGH A COMMON BLOCK TO SUBROUTINE CONVECT etc.
42  ! (Written by V.T.J. Phillips, 20-30/Jan/99)
43  ! ===========================================================================
44
45  ! line 1:  a flag (0 or 1) to decide whether P(F) = 1 or the general P(F)
46  ! is to be
47  ! used, followed by SCUT, which is the cut-off value of F in CONVECT
48  ! line 2:  blank
49  ! line 3:  the coefficients for the linear combination of P(F)s to
50  ! make the general P(F)
51  ! line 4:  blank
52  ! line 5:  gammas, Fmax for the cosh^2 component of P(F)
53  ! line 6:  blank
54  ! line 7:  alphas for the 1st irrational P(F)
55  ! line 8:  blank
56  ! line 9:  betas  for the 2nd irrational P(F)
57
58
59  ! c$$$        open(57,file='parameter_mix.data')
60  ! c$$$
61  ! c$$$        read(57,*) iflag_mix, scut
62  ! c$$$        read(57,*)
63  ! c$$$        if(iflag_mix .gt. 0) then
64  ! c$$$              read(57,*) qqa1, qqa2
65  ! c$$$              read(57,*)
66  ! c$$$              read(57,*) gammas, Fmax
67  ! c$$$              read(57,*)
68  ! c$$$              read(57,*) alphas
69  ! c$$$         endif
70  ! c$$$         close(57)
71
72
73  IF (iflag_mix>0) THEN
74
75    ! --      Normalize Pdf weights
76
77    sumcoef = qqa1 + qqa2
78    qqa1 = qqa1/sumcoef
79    qqa2 = qqa2/sumcoef
80
81    qcoef1max = qcoef1(fmax)
82    qcoef2max = qcoef2(fmax)
83
84    sigma = 0.
85    aire = 0.0
86    pdf = 0.0
87    mu = 0.0
88    df = 0.0001
89
90    ! do ff = 0.0 + df, 1.0 - 2.*df, df
91    ff = df
92    DO WHILE (ff<=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, *) pdf, qmix(ff), aire, ff
97      ff = ff + df
98    END DO
99
100    ! do ff=0.0+df,1.0 - 2.*df,df
101    ff = df
102    DO WHILE (ff<=1.0-2.*df)
103      pdf = (qmix(ff+df)-qmix(ff))*(1.-ff)/df
104      sigma = sigma + pdf*(ff-mu)*(ff-mu)*df
105      ff = ff + df
106    END DO
107    sigma = sqrt(sigma)
108
109    IF (abs(aire-1.0)>0.02) THEN
110      WRITE (lunout, *) 'WARNING:: AREA OF MIXING PDF IS::', aire
111      abort_message = ''
112      CALL abort_gcm(modname, abort_message, 1)
113    ELSE
114      PRINT *, 'Area, mean & std deviation are ::', aire, mu, sigma
115    END IF
116  END IF !  (iflag_mix .gt. 0)
117
118  RETURN
119END SUBROUTINE cv3_inip
Note: See TracBrowser for help on using the repository browser.