My Project
 All Classes Files Functions Variables Macros
nuage.F
Go to the documentation of this file.
1 ! $Id: nuage.F 1279 2009-12-10 09:02:56Z fairhead $
2 !
3  SUBROUTINE nuage (paprs, pplay,
4  . t, pqlwp, pclc, pcltau, pclemi,
5  . pch, pcl, pcm, pct, pctlwp,
6  e ok_aie,
7  e mass_solu_aero, mass_solu_aero_pi,
8  e bl95_b0, bl95_b1,
9  s cldtaupi, re, fl)
10  USE dimphy
11  IMPLICIT none
12 c======================================================================
13 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
14 c Objet: Calculer epaisseur optique et emmissivite des nuages
15 c======================================================================
16 c Arguments:
17 c t-------input-R-temperature
18 c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
19 c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
20 c ok_aie--input-L-apply aerosol indirect effect or not
21 c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
22 c mass_solu_aero_pi--input-R-dito, pre-industrial value
23 c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
24 c bl95_b1-input-R-a parameter, may be varied for tests ( -"- )
25 c
26 c cldtaupi-output-R-pre-industrial value of cloud optical thickness,
27 c needed for the diagnostics of the aerosol indirect
28 c radiative forcing (see radlwsw)
29 c re------output-R-Cloud droplet effective radius multiplied by fl [um]
30 c fl------output-R-Denominator to re, introduced to avoid problems in
31 c the averaging of the output. fl is the fraction of liquid
32 c water clouds within a grid cell
33 c
34 c pcltau--output-R-epaisseur optique des nuages
35 c pclemi--output-R-emissivite des nuages (0 a 1)
36 c======================================================================
37 C
38 #include "YOMCST.h"
39 c
40 cym#include "dimensions.h"
41 cym#include "dimphy.h"
42  REAL paprs(klon,klev+1), pplay(klon,klev)
43  REAL t(klon,klev)
44 c
45  REAL pclc(klon,klev)
46  REAL pqlwp(klon,klev)
47  REAL pcltau(klon,klev), pclemi(klon,klev)
48 c
49  REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
50 c
51  LOGICAL lo
52 c
53  REAL cetahb, cetamb
54  parameter(cetahb = 0.45, cetamb = 0.80)
55 C
56  INTEGER i, k
57  REAL zflwp, zradef, zfice, zmsac
58 c
59  REAL radius, rad_froid, rad_chaud, rad_chau1, rad_chau2
60  parameter(rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
61 ccc PARAMETER (rad_chaud=15.0, rad_froid=35.0)
62 c sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0)
63  REAL coef, coef_froi, coef_chau
64  parameter(coef_chau=0.13, coef_froi=0.09)
65  REAL seuil_neb, t_glace
66  parameter(seuil_neb=0.001, t_glace=273.0-15.0)
67  INTEGER nexpo ! exponentiel pour glace/eau
68  parameter(nexpo=6)
69 
70 cjq for the aerosol indirect effect
71 cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
72 cjq
73  LOGICAL ok_aie ! Apply AIE or not?
74 
75  REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols[ug m-3]
76  REAL mass_solu_aero_pi(klon, klev) ! - " - pre-industrial value
77  REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
78  REAL re(klon, klev) ! cloud droplet effective radius [um]
79  REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
80  REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
81 
82  REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
83 
84  REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
85 
86  REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
87 cjq-end
88 
89 ccc PARAMETER (nexpo=1)
90 c
91 c Calculer l'epaisseur optique et l'emmissivite des nuages
92 c
93  DO k = 1, klev
94  DO i = 1, klon
95  rad_chaud = rad_chau1
96  IF (k.LE.3) rad_chaud = rad_chau2
97 
98  pclc(i,k) = max(pclc(i,k), seuil_neb)
99  zflwp = 1000.*pqlwp(i,k)/rg/pclc(i,k)
100  . *(paprs(i,k)-paprs(i,k+1))
101  zfice = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
102  zfice = min(max(zfice,0.0),1.0)
103  zfice = zfice**nexpo
104 
105  IF (ok_aie) THEN
106  ! Formula "D" of Boucher and Lohmann, Tellus, 1995
107  !
108  cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
109  . log(max(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
110  ! Cloud droplet number concentration (CDNC) is restricted
111  ! to be within [20, 1000 cm^3]
112  !
113  cdnc(i,k)=min(1000.e6,max(20.e6,cdnc(i,k)))
114  cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
115  . log(max(mass_solu_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
116  cdnc_pi(i,k)=min(1000.e6,max(20.e6,cdnc_pi(i,k)))
117  !
118  !
119  ! air density: pplay(i,k) / (RD * zT(i,k))
120  ! factor 1.1: derive effective radius from volume-mean radius
121  ! factor 1000 is the water density
122  ! _chaud means that this is the CDR for liquid water clouds
123  !
124  rad_chaud =
125  . 1.1 * ( (pqlwp(i,k) * pplay(i,k) / (rd * t(i,k)) )
126  . / (4./3. * rpi * 1000. * cdnc(i,k)) )**(1./3.)
127  !
128  ! Convert to um. CDR shall be at least 3 um.
129  !
130  rad_chaud = max(rad_chaud*1.e6, 3.)
131 
132  ! For output diagnostics
133  !
134  ! Cloud droplet effective radius [um]
135  !
136  ! we multiply here with f * xl (fraction of liquid water
137  ! clouds in the grid cell) to avoid problems in the
138  ! averaging of the output.
139  ! In the output of IOIPSL, derive the real cloud droplet
140  ! effective radius as re/fl
141  !
142  fl(i,k) = pclc(i,k)*(1.-zfice)
143  re(i,k) = rad_chaud*fl(i,k)
144 
145  ! Pre-industrial cloud opt thickness
146  !
147  ! "radius" is calculated as rad_chaud above (plus the
148  ! ice cloud contribution) but using cdnc_pi instead of
149  ! cdnc.
150  radius = max(1.1e6 * ( (pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))
151  . / (4./3.*rpi*1000.*cdnc_pi(i,k)) )**(1./3.),
152  . 3.) * (1.-zfice) + rad_froid * zfice
153  cldtaupi(i,k) = 3.0/2.0 * zflwp / radius
154  .
155  ENDIF ! ok_aie
156 
157  radius = rad_chaud * (1.-zfice) + rad_froid * zfice
158  coef = coef_chau * (1.-zfice) + coef_froi * zfice
159  pcltau(i,k) = 3.0/2.0 * zflwp / radius
160  pclemi(i,k) = 1.0 - exp( - coef * zflwp)
161  lo = (pclc(i,k) .LE. seuil_neb)
162  IF (lo) pclc(i,k) = 0.0
163  IF (lo) pcltau(i,k) = 0.0
164  IF (lo) pclemi(i,k) = 0.0
165 
166  IF (.NOT.ok_aie) cldtaupi(i,k)=pcltau(i,k)
167  ENDDO
168  ENDDO
169 ccc DO k = 1, klev
170 ccc DO i = 1, klon
171 ccc t(i,k) = t(i,k)
172 ccc pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
173 ccc lo = pclc(i,k) .GT. (2.*1.e-5)
174 ccc zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
175 ccc . /(rg*pclc(i,k))
176 ccc zradef = 10.0 + (1.-sigs(k))*45.0
177 ccc pcltau(i,k) = 1.5 * zflwp / zradef
178 ccc zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
179 ccc zmsac = 0.13*(1.0-zfice) + 0.08*zfice
180 ccc pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
181 ccc if (.NOT.lo) pclc(i,k) = 0.0
182 ccc if (.NOT.lo) pcltau(i,k) = 0.0
183 ccc if (.NOT.lo) pclemi(i,k) = 0.0
184 ccc ENDDO
185 ccc ENDDO
186 cccccc print*, 'pas de nuage dans le rayonnement'
187 cccccc DO k = 1, klev
188 cccccc DO i = 1, klon
189 cccccc pclc(i,k) = 0.0
190 cccccc pcltau(i,k) = 0.0
191 cccccc pclemi(i,k) = 0.0
192 cccccc ENDDO
193 cccccc ENDDO
194 C
195 C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
196 C
197  DO i = 1, klon
198  pct(i)=1.0
199  pch(i)=1.0
200  pcm(i) = 1.0
201  pcl(i) = 1.0
202  pctlwp(i) = 0.0
203  ENDDO
204 C
205  DO k = klev, 1, -1
206  DO i = 1, klon
207  pctlwp(i) = pctlwp(i)
208  . + pqlwp(i,k)*(paprs(i,k)-paprs(i,k+1))/rg
209  pct(i) = pct(i)*(1.0-pclc(i,k))
210  if (pplay(i,k).LE.cetahb*paprs(i,1))
211  . pch(i) = pch(i)*(1.0-pclc(i,k))
212  if (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
213  . pplay(i,k).LE.cetamb*paprs(i,1))
214  . pcm(i) = pcm(i)*(1.0-pclc(i,k))
215  if (pplay(i,k).GT.cetamb*paprs(i,1))
216  . pcl(i) = pcl(i)*(1.0-pclc(i,k))
217  ENDDO
218  ENDDO
219 C
220  DO i = 1, klon
221  pct(i)=1.-pct(i)
222  pch(i)=1.-pch(i)
223  pcm(i)=1.-pcm(i)
224  pcl(i)=1.-pcl(i)
225  ENDDO
226 C
227  RETURN
228  END
229  SUBROUTINE diagcld1(paprs,pplay,rain,snow,kbot,ktop,
230  . diafra,dialiq)
231  use dimphy
232  IMPLICIT none
233 c
234 c Laurent Li (LMD/CNRS), le 12 octobre 1998
235 c (adaptation du code ECMWF)
236 c
237 c Dans certains cas, le schema pronostique des nuages n'est
238 c pas suffisament performant. On a donc besoin de diagnostiquer
239 c ces nuages. Je dois avouer que c'est une frustration.
240 c
241 cym#include "dimensions.h"
242 cym#include "dimphy.h"
243 #include "YOMCST.h"
244 c
245 c Arguments d'entree:
246  REAL paprs(klon,klev+1) ! pression (Pa) a inter-couche
247  REAL pplay(klon,klev) ! pression (Pa) au milieu de couche
248  REAL t(klon,klev) ! temperature (K)
249  REAL q(klon,klev) ! humidite specifique (Kg/Kg)
250  REAL rain(klon) ! pluie convective (kg/m2/s)
251  REAL snow(klon) ! neige convective (kg/m2/s)
252  INTEGER ktop(klon) ! sommet de la convection
253  INTEGER kbot(klon) ! bas de la convection
254 c
255 c Arguments de sortie:
256  REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
257  REAL dialiq(klon,klev) ! eau liquide nuageuse
258 c
259 c Constantes ajustables:
260  REAL canva, canvb, canvh
261  parameter(canva=2.0, canvb=0.3, canvh=0.4)
262  REAL cca, ccb, ccc
263  parameter(cca=0.125, ccb=1.5, ccc=0.8)
264  REAL ccfct, ccscal
265  parameter(ccfct=0.400)
266  parameter(ccscal=1.0e+11)
267  REAL cetahb, cetamb
268  parameter(cetahb=0.45, cetamb=0.80)
269  REAL cclwmr
270  parameter(cclwmr=1.e-04)
271  REAL zepscr
272  parameter(zepscr=1.0e-10)
273 c
274 c Variables locales:
275  INTEGER i, k
276  REAL zcc(klon)
277 c
278 c Initialisation:
279 c
280  DO k = 1, klev
281  DO i = 1, klon
282  diafra(i,k) = 0.0
283  dialiq(i,k) = 0.0
284  ENDDO
285  ENDDO
286 c
287  DO i = 1, klon ! Calculer la fraction nuageuse
288  zcc(i) = 0.0
289  IF((rain(i)+snow(i)).GT.0.) THEN
290  zcc(i)= cca * log(max(zepscr,(rain(i)+snow(i))*ccscal))-ccb
291  zcc(i)= min(ccc,max(0.0,zcc(i)))
292  ENDIF
293  ENDDO
294 c
295  DO i = 1, klon ! pour traiter les enclumes
296  diafra(i,ktop(i)) = max(diafra(i,ktop(i)),zcc(i)*ccfct)
297  IF ((zcc(i).GE.canvh) .AND.
298  . (pplay(i,ktop(i)).LE.cetahb*paprs(i,1)))
299  . diafra(i,ktop(i)) = max(diafra(i,ktop(i)),
300  . max(zcc(i)*ccfct,canva*(zcc(i)-canvb)))
301  dialiq(i,ktop(i))=cclwmr*diafra(i,ktop(i))
302  ENDDO
303 c
304  DO k = 1, klev ! nuages convectifs (sauf enclumes)
305  DO i = 1, klon
306  IF (k.LT.ktop(i) .AND. k.GE.kbot(i)) THEN
307  diafra(i,k)=max(diafra(i,k),zcc(i)*ccfct)
308  dialiq(i,k)=cclwmr*diafra(i,k)
309  ENDIF
310  ENDDO
311  ENDDO
312 c
313  RETURN
314  END
315  SUBROUTINE diagcld2(paprs,pplay,t,q, diafra,dialiq)
316  use dimphy
317  IMPLICIT none
318 c
319 cym#include "dimensions.h"
320 cym#include "dimphy.h"
321 #include "YOMCST.h"
322 c
323 c Arguments d'entree:
324  REAL paprs(klon,klev+1) ! pression (Pa) a inter-couche
325  REAL pplay(klon,klev) ! pression (Pa) au milieu de couche
326  REAL t(klon,klev) ! temperature (K)
327  REAL q(klon,klev) ! humidite specifique (Kg/Kg)
328 c
329 c Arguments de sortie:
330  REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
331  REAL dialiq(klon,klev) ! eau liquide nuageuse
332 c
333  REAL cetamb
334  parameter(cetamb=0.80)
335  REAL cloia, cloib, cloic, cloid
336  parameter(cloia=1.0e+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
337 ccc PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
338  REAL rgammas
339  parameter(rgammas=0.05)
340  REAL crhl
341  parameter(crhl=0.15)
342 ccc PARAMETER (CRHL=0.70)
343  REAL t_coup
344  parameter(t_coup=234.0)
345 c
346 c Variables locales:
347  INTEGER i, k, kb, invb(klon)
348  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
349  REAL zdelta, zcor
350 c
351 c Fonctions thermodynamiques:
352 #include "YOETHF.h"
353 #include "FCTTRE.h"
354 c
355 c Initialisation:
356 c
357  DO k = 1, klev
358  DO i = 1, klon
359  diafra(i,k) = 0.0
360  dialiq(i,k) = 0.0
361  ENDDO
362  ENDDO
363 c
364  DO i = 1, klon
365  invb(i) = klev
366  zdthmin(i)=0.0
367  ENDDO
368 
369  DO k = 2, klev/2-1
370  DO i = 1, klon
371  zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
372  . - rd * 0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i,k+1)
373  zdthdp = zdthdp * cloia
374  IF (pplay(i,k).GT.cetamb*paprs(i,1) .AND.
375  . zdthdp.LT.zdthmin(i) ) THEN
376  zdthmin(i) = zdthdp
377  invb(i) = k
378  ENDIF
379  ENDDO
380  ENDDO
381 
382  DO i = 1, klon
383  kb=invb(i)
384  IF (thermcep) THEN
385  zdelta=max(0.,sign(1.,rtt-t(i,kb)))
386  zqs= r2es*foeew(t(i,kb),zdelta)/pplay(i,kb)
387  zqs=min(0.5,zqs)
388  zcor=1./(1.-retv*zqs)
389  zqs=zqs*zcor
390  ELSE
391  IF (t(i,kb) .LT. t_coup) THEN
392  zqs = qsats(t(i,kb)) / pplay(i,kb)
393  ELSE
394  zqs = qsatl(t(i,kb)) / pplay(i,kb)
395  ENDIF
396  ENDIF
397  zcll = cloib * zdthmin(i) + cloic
398  zcll = min(1.0,max(0.0,zcll))
399  zrhb= q(i,kb)/zqs
400  IF (zcll.GT.0.0.AND.zrhb.LT.crhl)
401  . zcll=zcll*(1.-(crhl-zrhb)*cloid)
402  zcll=min(1.0,max(0.0,zcll))
403  diafra(i,kb) = max(diafra(i,kb),zcll)
404  dialiq(i,kb)= diafra(i,kb) * rgammas*zqs
405  ENDDO
406 c
407  RETURN
408  END