My Project
 All Classes Files Functions Variables Macros
newmicro.F
Go to the documentation of this file.
1 ! $Id: newmicro.F 1723 2013-02-04 13:25:36Z idelkadi $
2 
3 
4 !
5  SUBROUTINE newmicro (ok_cdnc, bl95_b0, bl95_b1,
6  . paprs, pplay,
7  . t, pqlwp, pclc, pcltau, pclemi,
8  . pch, pcl, pcm, pct, pctlwp,
9  . xflwp, xfiwp, xflwc, xfiwc,
10  . mass_solu_aero, mass_solu_aero_pi,
11  . pcldtaupi, re, fl, reliq, reice)
12 c
13  USE dimphy
14  USE phys_local_var_mod, only: scdnc,cldncl,reffclwtop,lcc,
15  . reffclws,reffclwc,cldnvi,lcc3d,
16  . lcc3dcon,lcc3dstra
17  USE phys_state_var_mod, only: rnebcon,clwcon
18  IMPLICIT none
19 c======================================================================
20 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
21 c O. Boucher (LMD/CNRS) mise a jour en 201212
22 c Objet: Calculer epaisseur optique et emmissivite des nuages
23 c======================================================================
24 c Arguments:
25 c ok_cdnc-input-L-flag pour calculer les rayons a partir des aerosols
26 c
27 c t-------input-R-temperature
28 c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere dans la partie nuageuse (kg/kg)
29 c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
30 c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
31 c mass_solu_aero_pi--input-R-ditto, pre-industrial value
32 c
33 c bl95_b0-input-R-a PARAMETER, may be varied for tests (s-sea, l-land)
34 c bl95_b1-input-R-a PARAMETER, may be varied for tests ( -"- )
35 c
36 c re------output-R-Cloud droplet effective radius multiplied by fl [um]
37 c fl------output-R-Denominator to re, introduced to avoid problems in
38 c the averaging of the output. fl is the fraction of liquid
39 c water clouds within a grid cell
40 c
41 c pcltau--output-R-epaisseur optique des nuages
42 c pclemi--output-R-emissivite des nuages (0 a 1)
43 c pcldtaupi-output-R-pre-industrial value of cloud optical thickness,
44 c
45 c pcl-output-R-2D low-level cloud cover
46 c pcm-output-R-2D mid-level cloud cover
47 c pch-output-R-2D high-level cloud cover
48 c pct-output-R-2D total cloud cover
49 c======================================================================
50 C
51 #include "YOMCST.h"
52 #include "nuage.h"
53 #include "radepsi.h"
54 #include "radopt.h"
55 
56 c choix de l'hypothese de recouvrememnt nuageuse
57  LOGICAL random, maximum_random, maximum
58  parameter(random=.false., maximum_random=.true., maximum=.false.)
59 c
60  LOGICAL, SAVE :: first=.true.
61 !$OMP THREADPRIVATE(FIRST)
62  INTEGER flag_max
63 c
64 c threshold PARAMETERs
65  REAL thres_tau,thres_neb
66  parameter(thres_tau=0.3, thres_neb=0.001)
67 c
68  REAL phase3d(klon, klev)
69  REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
70 c
71  REAL paprs(klon,klev+1)
72  REAL pplay(klon,klev)
73  REAL t(klon,klev)
74  REAL pclc(klon,klev)
75  REAL pqlwp(klon,klev)
76  REAL pcltau(klon,klev)
77  REAL pclemi(klon,klev)
78  REAL pcldtaupi(klon, klev)
79 c
80  REAL pct(klon)
81  REAL pcl(klon)
82  REAL pcm(klon)
83  REAL pch(klon)
84  REAL pctlwp(klon)
85 c
86  LOGICAL lo
87 c
88 !!Abderr modif JL mail du 19.01.2011 18:31
89 ! REAL cetahb, cetamb
90 ! PARAMETER (cetahb = 0.45, cetamb = 0.80)
91 ! Remplacer
92 ! cetahb*paprs(i,1) par prmhc
93 ! cetamb*paprs(i,1) par prlmc
94  REAL prmhc ! Pressure between medium and high level cloud in Pa
95  REAL prlmc ! Pressure between low and medium level cloud in Pa
96  parameter(prmhc = 440.*100., prlmc = 680.*100.)
97 C
98  INTEGER i, k
99  REAL xflwp(klon), xfiwp(klon)
100  REAL xflwc(klon,klev), xfiwc(klon,klev)
101 c
102  REAL radius
103 c
104  REAL coef_froi, coef_chau
105  parameter(coef_chau=0.13, coef_froi=0.09)
106 c
107  REAL seuil_neb
108  parameter(seuil_neb=0.001)
109 c
110  INTEGER nexpo ! exponentiel pour glace/eau
111  parameter(nexpo=6)
112 c PARAMETER (nexpo=1)
113 
114  REAL rel, tc, rei
115  REAL k_ice0, k_ice, df
116  parameter(k_ice0=0.005) ! units=m2/g
117  parameter(df=1.66) ! diffusivity factor
118 c
119 cjq for the aerosol indirect effect
120 cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
121 cjq
122  REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols [ug m-3]
123  REAL mass_solu_aero_pi(klon, klev) ! - " - (pre-industrial value)
124  REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
125  REAL re(klon, klev) ! cloud droplet effective radius [um]
126  REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
127  REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
128 
129  REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
130 
131  LOGICAL ok_cdnc
132  REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
133 
134 cjq-end
135 cIM cf. CR:parametres supplementaires
136  REAL zclear(klon)
137  REAL zcloud(klon)
138  REAL zcloudh(klon)
139  REAL zcloudm(klon)
140  REAL zcloudl(klon)
141  REAL rhodz(klon, klev) !--rho*dz pour la couche
142  REAL zrho(klon, klev) !--rho pour la couche
143  REAL dh(klon, klev) !--dz pour la couche
144  REAL zfice(klon, klev)
145  REAL rad_chaud(klon, klev) !--rayon pour les nuages chauds
146  REAL zflwp_var, zfiwp_var
147  REAL d_rei_dt
148 
149 ! Abderrahmane oct 2009
150  Real reliq(klon, klev), reice(klon, klev)
151 
152 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
153 ! FH : 2011/05/24
154 !
155 ! rei = ( rei_max - rei_min ) * T(°C) / 81.4 + rei_max
156 ! to be used for a temperature in celcius T(°C) < 0
157 ! rei=rei_min for T(°C) < -81.4
158 !
159 ! Calcul de la pente de la relation entre rayon effective des cristaux
160 ! et la température.
161 ! Pour retrouver les résultats numériques de la version d'origine,
162 ! on impose 0.71 quand on est proche de 0.71
163 c
164  d_rei_dt=(rei_max-rei_min)/81.4
165  if (abs(d_rei_dt-0.71)<1.e-4) d_rei_dt=0.71
166 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
167 c
168 c Calculer l'epaisseur optique et l'emmissivite des nuages
169 c IM inversion des DO
170 c
171  xflwp = 0.d0
172  xfiwp = 0.d0
173  xflwc = 0.d0
174  xfiwc = 0.d0
175 c
176  reliq=0.
177  reice=0.
178 c
179  DO k = 1, klev
180  DO i = 1, klon
181 c-layer calculation
182  rhodz(i,k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
183  zrho(i,k)=pplay(i,k)/t(i,k)/rd ! kg/m3
184  dh(i,k)=rhodz(i,k)/zrho(i,k) ! m
185 c-Fraction of ice in cloud using a linear transition
186  zfice(i,k) = 1.0 - (t(i,k)-t_glace_min) /
187  & (t_glace_max-t_glace_min)
188  zfice(i,k) = min(max(zfice(i,k),0.0),1.0)
189 c-IM Total Liquid/Ice water content
190  xflwc(i,k) = (1.-zfice(i,k))*pqlwp(i,k)
191  xfiwc(i,k) = zfice(i,k)*pqlwp(i,k)
192  ENDDO
193  ENDDO
194 
195  IF (ok_cdnc) THEN
196 c
197 c--we compute cloud properties as a function of the aerosol load
198 c
199  DO k = 1, klev
200  DO i = 1, klon
201 c
202 c Formula "D" of Boucher and Lohmann, Tellus, 1995
203 c Cloud droplet number concentration (CDNC) is restricted
204 c to be within [20, 1000 cm^3]
205 c
206 c--present-day case
207  cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
208  & log(max(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
209  cdnc(i,k)=min(1000.e6,max(20.e6,cdnc(i,k)))
210 c
211 c--pre-industrial case
212  cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
213  & log(max(mass_solu_aero_pi(i,k),1.e-4))/log(10.))
214  & *1.e6 !-m-3
215  cdnc_pi(i,k)=min(1000.e6,max(20.e6,cdnc_pi(i,k)))
216 c
217 c--present-day case
218  rad_chaud(i,k) =
219  & 1.1*((pqlwp(i,k)*pplay(i,k)/(rd * t(i,k)))
220  & /(4./3*rpi*1000.*cdnc(i,k)) )**(1./3.)
221  rad_chaud(i,k) = max(rad_chaud(i,k) * 1.e6, 5.)
222 c
223 c--pre-industrial case
224  radius =
225  & 1.1*((pqlwp(i,k)*pplay(i,k)/(rd * t(i,k)))
226  & /(4./3.*rpi*1000.*cdnc_pi(i,k)))**(1./3.)
227  radius = max(radius*1.e6, 5.)
228 c
229 c--pre-industrial case
230 c--liquid/ice cloud water paths:
231  IF (pclc(i,k) .LE. seuil_neb) THEN
232 c
233  pcldtaupi(i,k) = 0.0
234 c
235  ELSE
236 c
237  zflwp_var= 1000.*(1.-zfice(i,k))*pqlwp(i,k)/pclc(i,k)
238  & *rhodz(i,k)
239  zfiwp_var= 1000.*zfice(i,k)*pqlwp(i,k)/pclc(i,k)
240  & *rhodz(i,k)
241  tc = t(i,k)-273.15
242  rei = d_rei_dt*tc + rei_max
243  if (tc.le.-81.4) rei = rei_min
244 c
245 c-- cloud optical thickness :
246 c [for liquid clouds, traditional formula,
247 c for ice clouds, Ebert & Curry (1992)]
248 c
249  if (zflwp_var.eq.0.) radius = 1.
250  if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
251  pcldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius
252  & + zfiwp_var * (3.448e-03 + 2.431/rei)
253 c
254  ENDIF
255 c
256  ENDDO
257  ENDDO
258 c
259  ELSE !--not ok_cdnc
260 c
261 c-prescribed cloud droplet radius
262 c
263  DO k = 1, min(3,klev)
264  DO i = 1, klon
265  rad_chaud(i,k) = rad_chau2
266  ENDDO
267  ENDDO
268  DO k = min(3,klev)+1, klev
269  DO i = 1, klon
270  rad_chaud(i,k) = rad_chau1
271  ENDDO
272  ENDDO
273 
274  ENDIF !--ok_cdnc
275 c
276 c--computation of cloud optical depth and emissivity
277 c--in the general case
278 c
279  DO k = 1, klev
280  DO i = 1, klon
281 c
282  IF (pclc(i,k) .LE. seuil_neb) THEN
283 c
284 c effective cloud droplet radius (microns) for liquid water clouds:
285 c For output diagnostics cloud droplet effective radius [um]
286 c we multiply here with f * xl (fraction of liquid water
287 c clouds in the grid cell) to avoid problems in the averaging of the output.
288 c In the output of IOIPSL, derive the REAL cloud droplet
289 c effective radius as re/fl
290 c
291  fl(i,k) = seuil_neb*(1.-zfice(i,k))
292  re(i,k) = rad_chaud(i,k)*fl(i,k)
293  rel = 0.
294  rei = 0.
295  pclc(i,k) = 0.0
296  pcltau(i,k) = 0.0
297  pclemi(i,k) = 0.0
298 c
299  ELSE
300 c
301 c -- liquid/ice cloud water paths:
302 
303  zflwp_var= 1000.*(1.-zfice(i,k))*pqlwp(i,k)/pclc(i,k)
304  & *rhodz(i,k)
305  zfiwp_var= 1000.*zfice(i,k)*pqlwp(i,k)/pclc(i,k)
306  & *rhodz(i,k)
307 c
308 c effective cloud droplet radius (microns) for liquid water clouds:
309 c For output diagnostics cloud droplet effective radius [um]
310 c we multiply here with f * xl (fraction of liquid water
311 c clouds in the grid cell) to avoid problems in the averaging of the output.
312 c In the output of IOIPSL, derive the REAL cloud droplet
313 c effective radius as re/fl
314 c
315  fl(i,k) = pclc(i,k)*(1.-zfice(i,k))
316  re(i,k) = rad_chaud(i,k)*fl(i,k)
317 c
318  rel = rad_chaud(i,k)
319 c
320 c for ice clouds: as a function of the ambiant temperature
321 c [formula used by Iacobellis and Somerville (2000), with an
322 c asymptotical value of 3.5 microns at T<-81.4 C added to be
323 c consistent with observations of Heymsfield et al. 1986]:
324 c 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as rei_max=61.29
325 c
326  tc = t(i,k)-273.15
327  rei = d_rei_dt*tc + rei_max
328  if (tc.le.-81.4) rei = rei_min
329 c
330 c-- cloud optical thickness :
331 c [for liquid clouds, traditional formula,
332 c for ice clouds, Ebert & Curry (1992)]
333 c
334  if (zflwp_var.eq.0.) rel = 1.
335  if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
336  pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel )
337  & + zfiwp_var * (3.448e-03 + 2.431/rei)
338 c
339 c -- cloud infrared emissivity:
340 c [the broadband infrared absorption coefficient is PARAMETERized
341 c as a function of the effective cld droplet radius]
342 c Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
343 c
344  k_ice = k_ice0 + 1.0/rei
345 c
346  pclemi(i,k) = 1.0
347  & - exp( -coef_chau*zflwp_var - df*k_ice*zfiwp_var)
348 c
349  ENDIF
350 c
351  reliq(i,k)=rel
352  reice(i,k)=rei
353 c
354  xflwp(i) = xflwp(i)+ xflwc(i,k) * rhodz(i,k)
355  xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * rhodz(i,k)
356 c
357  ENDDO
358  ENDDO
359 c
360 c--if cloud droplet radius is fixed, then pcldtaupi=pcltau
361 c
362  IF (.NOT.ok_cdnc) THEN
363  DO k = 1, klev
364  DO i = 1, klon
365  pcldtaupi(i,k)=pcltau(i,k)
366  ENDDO
367  ENDDO
368  ENDIF
369 C
370 C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
371 c IM cf. CR:test: calcul prenant ou non en compte le recouvrement
372 c initialisations
373 c
374  DO i=1,klon
375  zclear(i)=1.
376  zcloud(i)=0.
377  zcloudh(i)=0.
378  zcloudm(i)=0.
379  zcloudl(i)=0.
380  pch(i)=1.0
381  pcm(i) = 1.0
382  pcl(i) = 1.0
383  pctlwp(i) = 0.0
384  ENDDO
385 C
386 c--calculation of liquid water path
387 c
388  DO k = klev, 1, -1
389  DO i = 1, klon
390  pctlwp(i) = pctlwp(i)+ pqlwp(i,k)*rhodz(i,k)
391  ENDDO
392  ENDDO
393 c
394 c--calculation of cloud properties with cloud overlap
395 c
396  IF (novlp.EQ.1) THEN
397  DO k = klev, 1, -1
398  DO i = 1, klon
399  zclear(i)=zclear(i)*(1.-max(pclc(i,k),zcloud(i)))
400  & /(1.-min(REAL(zcloud(i), kind=8),1.-zepsec))
401  pct(i)=1.-zclear(i)
402  IF (paprs(i,k).LT.prmhc) THEN
403  pch(i) = pch(i)*(1.-max(pclc(i,k),zcloudh(i)))
404  & /(1.-min(REAL(zcloudh(i), kind=8),1.-zepsec))
405  zcloudh(i)=pclc(i,k)
406  ELSE IF (paprs(i,k).GE.prmhc .AND.
407  & paprs(i,k).LT.prlmc) THEN
408  pcm(i) = pcm(i)*(1.-max(pclc(i,k),zcloudm(i)))
409  & /(1.-min(REAL(zcloudm(i), kind=8),1.-zepsec))
410  zcloudm(i)=pclc(i,k)
411  ELSE IF (paprs(i,k).GE.prlmc) THEN
412  pcl(i) = pcl(i)*(1.-max(pclc(i,k),zcloudl(i)))
413  & /(1.-min(REAL(zcloudl(i), kind=8),1.-zepsec))
414  zcloudl(i)=pclc(i,k)
415  endif
416  zcloud(i)=pclc(i,k)
417  ENDDO
418  ENDDO
419  ELSE IF (novlp.EQ.2) THEN
420  DO k = klev, 1, -1
421  DO i = 1, klon
422  zcloud(i)=max(pclc(i,k),zcloud(i))
423  pct(i)=zcloud(i)
424  IF (paprs(i,k).LT.prmhc) THEN
425  pch(i) = min(pclc(i,k),pch(i))
426  ELSE IF (paprs(i,k).GE.prmhc .AND.
427  & paprs(i,k).LT.prlmc) THEN
428  pcm(i) = min(pclc(i,k),pcm(i))
429  ELSE IF (paprs(i,k).GE.prlmc) THEN
430  pcl(i) = min(pclc(i,k),pcl(i))
431  endif
432  ENDDO
433  ENDDO
434  ELSE IF (novlp.EQ.3) THEN
435  DO k = klev, 1, -1
436  DO i = 1, klon
437  zclear(i)=zclear(i)*(1.-pclc(i,k))
438  pct(i)=1-zclear(i)
439  IF (paprs(i,k).LT.prmhc) THEN
440  pch(i) = pch(i)*(1.0-pclc(i,k))
441  ELSE IF (paprs(i,k).GE.prmhc .AND.
442  & paprs(i,k).LT.prlmc) THEN
443  pcm(i) = pcm(i)*(1.0-pclc(i,k))
444  ELSE IF (paprs(i,k).GE.prlmc) THEN
445  pcl(i) = pcl(i)*(1.0-pclc(i,k))
446  endif
447  ENDDO
448  ENDDO
449  ENDIF
450 C
451  DO i = 1, klon
452  pch(i)=1.-pch(i)
453  pcm(i)=1.-pcm(i)
454  pcl(i)=1.-pcl(i)
455  ENDDO
456 c
457 c ========================================================
458 c DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
459 c ========================================================
460 c change by Nicolas Yan (LSCE)
461 c Cloud Droplet Number Concentration (CDNC) : 3D variable
462 c Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
463 c Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
464 c Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
465 c Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
466 c
467  IF (ok_cdnc) THEN
468 c
469  DO k = 1, klev
470  DO i = 1, klon
471  phase3d(i,k)=1-zfice(i,k)
472  IF (pclc(i,k) .LE. seuil_neb) THEN
473  lcc3d(i,k)=seuil_neb*phase3d(i,k)
474  ELSE
475  lcc3d(i,k)=pclc(i,k)*phase3d(i,k)
476  ENDIF
477  scdnc(i,k)=lcc3d(i,k)*cdnc(i,k) ! m-3
478  ENDDO
479  ENDDO
480 c
481  DO i=1,klon
482  lcc(i)=0.
483  reffclwtop(i)=0.
484  cldncl(i)=0.
485  IF(random .OR. maximum_random) tcc(i) = 1.
486  IF(maximum) tcc(i) = 0.
487  ENDDO
488 c
489  DO i=1,klon
490  DO k=klev-1,1,-1 !From TOA down
491 c
492  ! Test, if the cloud optical depth exceeds the necessary
493  ! threshold:
494 
495  IF (pcltau(i,k).GT.thres_tau
496  . .AND. pclc(i,k).GT.thres_neb) THEN
497 
498  IF (maximum) THEN
499  IF (first) THEN
500  write(*,*)'Hypothese de recouvrement: MAXIMUM'
501  first=.false.
502  ENDIF
503  flag_max= -1.
504  ftmp(i) = max(tcc(i),pclc(i,k))
505  ENDIF
506 
507  IF (random) THEN
508  IF (first) THEN
509  write(*,*)'Hypothese de recouvrement: RANDOM'
510  first=.false.
511  ENDIF
512  flag_max= 1.
513  ftmp(i) = tcc(i) * (1-pclc(i,k))
514  ENDIF
515 
516  IF (maximum_random) THEN
517  IF (first) THEN
518  write(*,*)
519 'Hypothese de recouvrement: MAXIMUM_ . RANDOM'
520  first=.false.
521  ENDIF
522  flag_max= 1.
523  ftmp(i) = tcc(i) *
524  . (1. - max(pclc(i,k),pclc(i,k+1))) /
525  . (1. - min(pclc(i,k+1),1.-thres_neb))
526  ENDIF
527 c Effective radius of cloud droplet at top of cloud (m)
528  reffclwtop(i) = reffclwtop(i) + rad_chaud(i,k) *
529  . 1.0e-06 * phase3d(i,k) * ( tcc(i) - ftmp(i))*flag_max
530 c CDNC at top of cloud (m-3)
531  cldncl(i) = cldncl(i) + cdnc(i,k) * phase3d(i,k) *
532  . (tcc(i) - ftmp(i))*flag_max
533 c Liquid Cloud Content at top of cloud
534  lcc(i) = lcc(i) + phase3d(i,k) * (tcc(i)-ftmp(i))*
535  . flag_max
536 c Total Cloud Content at top of cloud
537  tcc(i)=ftmp(i)
538 
539  ENDIF ! is there a visible, not-too-small cloud?
540  ENDDO ! loop over k
541 c
542  IF (random .OR. maximum_random) tcc(i)=1.-tcc(i)
543 c
544  ENDDO ! loop over i
545 
546 !! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC REFFCLWS)
547  DO i = 1, klon
548  DO k = 1, klev
549 ! Weight to be used for outputs: eau_liquide*couverture nuageuse
550  lcc3dcon(i,k) =rnebcon(i,k)*phase3d(i,k)*clwcon(i,k) ! eau liquide convective
551  lcc3dstra(i,k)=pclc(i,k)*pqlwp(i,k)*phase3d(i,k)
552  lcc3dstra(i,k)=lcc3dstra(i,k)-lcc3dcon(i,k) ! eau liquide stratiforme
553  lcc3dstra(i,k)=max(lcc3dstra(i,k),0.0)
554 ! Compute cloud droplet radius as above in meter
555  radius=1.1*((pqlwp(i,k)*pplay(i,k)/(rd * t(i,k)))
556  & /(4./3*rpi*1000.*cdnc(i,k)) )**(1./3.)
557  radius=max(radius, 5.e-6)
558 ! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
559  reffclwc(i,k)=radius
560  reffclwc(i,k)=reffclwc(i,k)*lcc3dcon(i,k)
561 ! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
562  reffclws(i,k)=radius
563  reffclws(i,k)=reffclws(i,k)*lcc3dstra(i,k)
564  ENDDO !klev
565  ENDDO !klon
566 c
567 c Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
568 c
569  DO i = 1, klon
570  cldnvi(i)=0.
571  lcc_integrat(i)=0.
572  height(i)=0.
573  DO k = 1, klev
574  cldnvi(i)=cldnvi(i)+cdnc(i,k)*lcc3d(i,k)*dh(i,k)
575  lcc_integrat(i)=lcc_integrat(i)+lcc3d(i,k)*dh(i,k)
576  height(i)=height(i)+dh(i,k)
577  ENDDO ! klev
578  lcc_integrat(i)=lcc_integrat(i)/height(i)
579  IF (lcc_integrat(i) .LE. 1.0e-03) THEN
580  cldnvi(i)=cldnvi(i)*lcc(i)/seuil_neb
581  ELSE
582  cldnvi(i)=cldnvi(i)*lcc(i)/lcc_integrat(i)
583  ENDIF
584  ENDDO ! klon
585 
586  DO i = 1, klon
587  DO k = 1, klev
588  IF (scdnc(i,k) .LE. 0.0) scdnc(i,k)=0.0
589  IF (reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
590  IF (reffclwc(i,k) .LE. 0.0) reffclwc(i,k)=0.0
591  IF (lcc3d(i,k) .LE. 0.0) lcc3d(i,k)=0.0
592  IF (lcc3dcon(i,k) .LE. 0.0) lcc3dcon(i,k)=0.0
593  IF (lcc3dstra(i,k) .LE. 0.0) lcc3dstra(i,k)=0.0
594  ENDDO
595  IF (reffclwtop(i) .LE. 0.0) reffclwtop(i)=0.0
596  IF (cldncl(i) .LE. 0.0) cldncl(i)=0.0
597  IF (cldnvi(i) .LE. 0.0) cldnvi(i)=0.0
598  IF (lcc(i) .LE. 0.0) lcc(i)=0.0
599  ENDDO
600 c
601  ENDIF !ok_cdnc
602 c
603  RETURN
604 c
605  END