source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/cloudth.F90 @ 3069

Last change on this file since 3069 was 2160, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand

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