source: LMDZ5/trunk/libf/phylmd/cloudth.F90 @ 1992

Last change on this file since 1992 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • 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.3 KB
RevLine 
[1399]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
[1411]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) 
[1399]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
[1411]205      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
[1399]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
Note: See TracBrowser for help on using the repository browser.