source: LMDZ5/trunk/libf/phylmd/cv3_inicp.F90 @ 2201

Last change on this file since 2201 was 2197, checked in by Ehouarn Millour, 10 years ago

Added 'implicit none' statements and proper variable definitions where they were missing.
EM

  • 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.4 KB
Line 
1SUBROUTINE 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_gcm(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
118END SUBROUTINE cv3_inicp
Note: See TracBrowser for help on using the repository browser.