My Project
 All Classes Files Functions Variables Macros
cloudth.F90
Go to the documentation of this file.
1  SUBROUTINE cloudth(ngrid,klev,ind2, &
2  & ztv,po,zqta,fraca, &
3  & qcloud,ctot,zpspsk,paprs,ztla,zthl, &
4  & ratqs,zqs,t)
5 
6 
7  IMPLICIT NONE
8 
9 
10 !===========================================================================
11 ! Auteur : Arnaud Octavio Jam (LMD/CNRS)
12 ! Date : 25 Mai 2010
13 ! Objet : calcule les valeurs de qc et rneb dans les thermiques
14 !===========================================================================
15 
16 
17 #include "YOMCST.h"
18 #include "YOETHF.h"
19 #include "FCTTRE.h"
20 #include "iniprint.h"
21 #include "thermcell.h"
22 
23  INTEGER itap,ind1,ind2
24  INTEGER ngrid,klev,klon,l,ig
25 
26  REAL ztv(ngrid,klev)
27  REAL po(ngrid)
28  REAL zqenv(ngrid)
29  REAL zqta(ngrid,klev)
30 
31  REAL fraca(ngrid,klev+1)
32  REAL zpspsk(ngrid,klev)
33  REAL paprs(ngrid,klev+1)
34  REAL ztla(ngrid,klev)
35  REAL zthl(ngrid,klev)
36 
37  REAL zqsatth(ngrid,klev)
38  REAL zqsatenv(ngrid,klev)
39 
40 
41  REAL sigma1(ngrid,klev)
42  REAL sigma2(ngrid,klev)
43  REAL qlth(ngrid,klev)
44  REAL qlenv(ngrid,klev)
45  REAL qltot(ngrid,klev)
46  REAL cth(ngrid,klev)
47  REAL cenv(ngrid,klev)
48  REAL ctot(ngrid,klev)
49  REAL rneb(ngrid,klev)
50  REAL t(ngrid,klev)
51  REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
52  REAL rdd,cppd,lv
53  REAL alth,alenv,ath,aenv
54  REAL sth,senv,sigma1s,sigma2s,xth,xenv
55  REAL tbef,zdelta,qsatbef,zcor
56  REAL alpha,qlbef
57  REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
58 
59  REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
60  REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
61  REAL zqs(ngrid), qcloud(ngrid)
62  REAL erf
63 
64 
65 
66 
67 
68 ! print*,ngrid,klev,ind1,ind2,ztv(ind1,ind2),po(ind1),zqta(ind1,ind2), &
69 ! & fraca(ind1,ind2),zpspsk(ind1,ind2),paprs(ind1,ind2),ztla(ind1,ind2),zthl(ind1,ind2), &
70 ! & 'verif'
71 
72 
73 ! LOGICAL active(ngrid)
74 
75 !-----------------------------------------------------------------------------------------------------------------
76 ! Initialisation des variables réelles
77 !-----------------------------------------------------------------------------------------------------------------
78  sigma1(:,:)=0.
79  sigma2(:,:)=0.
80  qlth(:,:)=0.
81  qlenv(:,:)=0.
82  qltot(:,:)=0.
83  rneb(:,:)=0.
84  qcloud(:)=0.
85  cth(:,:)=0.
86  cenv(:,:)=0.
87  ctot(:,:)=0.
88  qsatmmussig1=0.
89  qsatmmussig2=0.
90  rdd=287.04
91  cppd=1005.7
92  pi=3.14159
93  lv=2.5e6
94  sqrt2pi=sqrt(2.*pi)
95 
96 
97 
98 !------------------------------------------------------------------------------------------------------------------
99 ! Calcul de la fraction du thermique et des écart-types des distributions
100 !------------------------------------------------------------------------------------------------------------------
101  do ind1=1,ngrid
102 
103  if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
104 
105  zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
106 
107 
108 ! zqenv(ind1)=po(ind1)
109  tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
110  zdelta=max(0.,sign(1.,rtt-tbef))
111  qsatbef= r2es * foeew(tbef,zdelta)/paprs(ind1,ind2)
112  qsatbef=min(0.5,qsatbef)
113  zcor=1./(1.-retv*qsatbef)
114  qsatbef=qsatbef*zcor
115  zqsatenv(ind1,ind2)=qsatbef
116 
117 
118 
119 
120  alenv=(0.622*lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
121  aenv=1./(1.+(alenv*lv/cppd))
122  senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
123 
124 
125 
126 
127  tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
128  zdelta=max(0.,sign(1.,rtt-tbef))
129  qsatbef= r2es * foeew(tbef,zdelta)/paprs(ind1,ind2)
130  qsatbef=min(0.5,qsatbef)
131  zcor=1./(1.-retv*qsatbef)
132  qsatbef=qsatbef*zcor
133  zqsatth(ind1,ind2)=qsatbef
134 
135  alth=(0.622*lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)
136  ath=1./(1.+(alth*lv/cppd))
137  sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
138 
139 
140 
141 !-----------------------------------------------------------------------------------------------------------------
142 ! Calcul des écart-types pour s
143 !-----------------------------------------------------------------------------------------------------------------
144 
145  sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
146  sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
147 
148 
149 !-----------------------------------------------------------------------------------------------------------------
150 ! Calcul de l'eau condensée et de la couverture nuageuse
151 !-----------------------------------------------------------------------------------------------------------------
152  sqrt2pi=sqrt(2.*pi)
153  xth=sth/(sqrt(2.)*sigma2s)
154  xenv=senv/(sqrt(2.)*sigma1s)
155  cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
156  cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
157  ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
158 ! ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
159 
160 
161 
162  qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
163  qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
164  qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
165 ! qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
166 
167 
168 ! print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
169 
170 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
171  if (ctot(ind1,ind2).lt.1.e-10) then
172  ctot(ind1,ind2)=0.
173  qcloud(ind1)=zqsatenv(ind1,ind2)
174 
175  else
176 
177  ctot(ind1,ind2)=ctot(ind1,ind2)
178  qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
179 
180  endif
181 
182 
183 ! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
184 
185 
186  else ! gaussienne environnement seule
187 
188  zqenv(ind1)=po(ind1)
189  tbef=t(ind1,ind2)
190  zdelta=max(0.,sign(1.,rtt-tbef))
191  qsatbef= r2es * foeew(tbef,zdelta)/paprs(ind1,ind2)
192  qsatbef=min(0.5,qsatbef)
193  zcor=1./(1.-retv*qsatbef)
194  qsatbef=qsatbef*zcor
195  zqsatenv(ind1,ind2)=qsatbef
196 
197 
198 ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
199  zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
200  alenv=(0.622*lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
201  aenv=1./(1.+(alenv*lv/cppd))
202  senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
203 
204 
205  sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
206 
207  sqrt2pi=sqrt(2.*pi)
208  xenv=senv/(sqrt(2.)*sigma1s)
209  ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
210  qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
211 
212  if (ctot(ind1,ind2).lt.1.e-3) then
213  ctot(ind1,ind2)=0.
214  qcloud(ind1)=zqsatenv(ind1,ind2)
215 
216  else
217 
218  ctot(ind1,ind2)=ctot(ind1,ind2)
219  qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
220 
221  endif
222 
223 
224 
225 
226 
227 
228  endif
229  enddo
230 
231  return
232  end
233 
234 
235 
236 
237 
238 
239 
240 
241 
242 
243