source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cloudth.F90 @ 5186

Last change on this file since 5186 was 1399, checked in by idelkadi, 14 years ago

Modifications relatives a l'inclusion des nouveaux nuages d'Arnaud Jam

  • calltherm.F90 : remontee de variables suplementaires.

Passage de ztv, zpspsk, zthla, zthl en argument

  • fisrtilp.F : le programme de nuages avec appel a cloudth

Ajout des arguments : zqta, fraca, ztv, zpspsk, ztla, ztlh,
iflag_cldcon

  • physiq.F : appel modifie a calltherm et fisrtilp
  • thermcell_main.F90 : changement dans les flag_thermals_ed
  • thermcell_plume.F90 : version modifiee entrainement et detrainement

(on retrouve la precendente pour iflag_thermals_ed=10)

  • cloudth.F90 : PDF d'eau bi-gaussiennes
File size: 8.0 KB
Line 
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!      alpha=0.5*(fraca(ind1,ind2)+fraca(ind1,ind2-1))
145!      sigma1s=(Max(1.2*fraca(ind1,ind2)*(sth-senv)**2,(senv/100)**2))**0.5
146!      sigma2s=((0.021/(fraca(ind1,ind2)+0.0)**0.5)*(sth-senv)**2+(sth/100)**2)**0.5
147!      sigma1s=(1.5**0.5)*(fraca(ind1,ind2)**0.65)*(sth-senv)+0.000045
148!      sigma2s=0.1265*(sth-senv)/(fraca(ind1,ind2)+0.0)**0.35+0.000045     
149
150      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.00003
151      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
152
153!      sigma1s=(1.5**0.5)*(alpha**0.65)*(sth-senv)+0.000045
154!      sigma2s=0.1265*(sth-senv)/alpha**0.35+0.000045
155     
156!      sigma1s=(2.8**0.5)*(0.1**0.7)*(sth-senv)+0.00002
157!      sigma2s=((0.126/(0.1)**0.3)*(sth-senv))+0.00002
158
159
160 
161!-----------------------------------------------------------------------------------------------------------------
162! Calcul de l'eau condensée et de la couverture nuageuse
163!-----------------------------------------------------------------------------------------------------------------
164      sqrt2pi=sqrt(2.*pi)
165      xth=sth/(sqrt(2.)*sigma2s)
166      xenv=senv/(sqrt(2.)*sigma1s)
167      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
168      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
169      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
170!      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
171
172
173
174      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
175      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
176      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
177!      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
178     
179
180!      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
181
182!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
183      if (ctot(ind1,ind2).lt.1.e-10) then
184      ctot(ind1,ind2)=0.
185      qcloud(ind1)=zqsatenv(ind1,ind2)
186
187      else   
188               
189      ctot(ind1,ind2)=ctot(ind1,ind2)
190      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
191
192      endif                           
193     
194         
195!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
196
197
198      else  ! gaussienne environnement seule
199     
200      zqenv(ind1)=po(ind1)
201      Tbef=t(ind1,ind2)
202      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
203      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
204      qsatbef=MIN(0.5,qsatbef)
205      zcor=1./(1.-retv*qsatbef)
206      qsatbef=qsatbef*zcor
207      zqsatenv(ind1,ind2)=qsatbef
208     
209
210!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
211      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
212      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
213      aenv=1./(1.+(alenv*Lv/cppd))
214      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
215     
216!      if (zqenv(ind1).gt.0.005) then
217      sigma1s=0.005*zqenv(ind1)
218!      else
219!      sigma1s=0.0005
220!      endif
221
222!      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
223
224
225
226!      sigma1s=0.00003
227      sqrt2pi=sqrt(2.*pi)
228      xenv=senv/(sqrt(2.)*sigma1s)
229      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
230      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
231     
232      if (ctot(ind1,ind2).lt.1.e-3) then
233      ctot(ind1,ind2)=0.
234      qcloud(ind1)=zqsatenv(ind1,ind2)
235
236      else   
237               
238      ctot(ind1,ind2)=ctot(ind1,ind2)
239      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
240
241      endif   
242 
243 
244 
245 
246 
247 
248      endif   
249      enddo
250     
251      return
252      end
253
254
255
256
257
258                                                                           
259
260
261
262
263
Note: See TracBrowser for help on using the repository browser.