source: LMDZ5/trunk/libf/phylmd/newmicro.F @ 1712

Last change on this file since 1712 was 1712, checked in by idelkadi, 11 years ago

Corrections dans newmicro.F pour completer les modifications de O. Boucher sur ok_ade/ok_aie avec correction des diagnostiques CMIP5 :
Rajout d'un nouveau flag ok_cdnc (ok cloud droplet number concentration)
Dans le cas sans aérosols, nous avons flag_aerosol=0, ok_cdnc=n, ok_ade=n et ok_aie=n
Dans le cas avec aérosols, nous avons flag_aerosol=6, ok_cdnc=y, ok_ade=y et ok_aie=y

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.1 KB
Line 
1! $Id: newmicro.F 1712 2013-01-18 14:07:29Z 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)
12c
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
19c======================================================================
20c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
21c            O.   Boucher (LMD/CNRS) mise a jour en 201212
22c Objet: Calculer epaisseur optique et emmissivite des nuages
23c======================================================================
24c Arguments:
25c ok_cdnc-input-L-flag pour calculer les rayons a partir des aerosols
26c
27c t-------input-R-temperature
28c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere dans la partie nuageuse (kg/kg)
29c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
30c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
31c mass_solu_aero_pi--input-R-ditto, pre-industrial value
32c
33c bl95_b0-input-R-a PARAMETER, may be varied for tests (s-sea, l-land)
34c bl95_b1-input-R-a PARAMETER, may be varied for tests (    -"-      )
35c     
36c re------output-R-Cloud droplet effective radius multiplied by fl [um]
37c fl------output-R-Denominator to re, introduced to avoid problems in
38c                  the averaging of the output. fl is the fraction of liquid
39c                  water clouds within a grid cell           
40c
41c pcltau--output-R-epaisseur optique des nuages
42c pclemi--output-R-emissivite des nuages (0 a 1)
43c pcldtaupi-output-R-pre-industrial value of cloud optical thickness,
44c
45c pcl-output-R-2D low-level cloud cover
46c pcm-output-R-2D mid-level cloud cover
47c pch-output-R-2D high-level cloud cover
48c pct-output-R-2D total cloud cover
49c======================================================================
50C
51#include "YOMCST.h"
52#include "nuage.h"
53#include "radepsi.h"
54#include "radopt.h"
55
56c choix de l'hypothese de recouvrememnt nuageuse
57      LOGICAL RANDOM, MAXIMUM_RANDOM, MAXIMUM
58      PARAMETER (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.)
59c
60      LOGICAL, SAVE :: FIRST=.TRUE.
61!$OMP THREADPRIVATE(FIRST)
62      INTEGER flag_max
63c
64c threshold PARAMETERs
65      REAL thres_tau,thres_neb
66      PARAMETER (thres_tau=0.3, thres_neb=0.001)
67c
68      REAL phase3d(klon, klev)
69      REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
70c
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)
79c
80      REAL pct(klon)
81      REAL pcl(klon)
82      REAL pcm(klon)
83      REAL pch(klon)
84      REAL pctlwp(klon)
85c
86      LOGICAL lo
87c
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.)
97C
98      INTEGER i, k
99      REAL xflwp(klon), xfiwp(klon)
100      REAL xflwc(klon,klev), xfiwc(klon,klev)
101c
102      REAL radius
103c
104      REAL coef_froi, coef_chau
105      PARAMETER (coef_chau=0.13, coef_froi=0.09)
106c
107      REAL seuil_neb
108      PARAMETER (seuil_neb=0.001)
109c
110      INTEGER nexpo ! exponentiel pour glace/eau
111      PARAMETER (nexpo=6)
112c      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
118c
119cjq for the aerosol indirect effect
120cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
121cjq     
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     
134cjq-end   
135cIM 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
163c
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!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
167c
168c Calculer l'epaisseur optique et l'emmissivite des nuages
169c     IM inversion des DO
170c
171      xflwp = 0.d0
172      xfiwp = 0.d0
173      xflwc = 0.d0
174      xfiwc = 0.d0
175c
176      reliq=0.
177      reice=0.
178c
179      DO k = 1, klev
180        DO i = 1, klon
181c-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
185c-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)
189c-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
196c
197c--we compute cloud properties as a function of the aerosol load
198c
199        DO k = 1, klev
200            DO i = 1, klon
201c
202c Formula "D" of Boucher and Lohmann, Tellus, 1995
203c Cloud droplet number concentration (CDNC) is restricted
204c to be within [20, 1000 cm^3]
205c
206c--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)))
210c
211c--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)))
216c
217c--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.)
222c
223c--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.)
228c
229c--pre-industrial case
230c--liquid/ice cloud water paths:
231                IF (pclc(i,k) .LE. seuil_neb) THEN
232c
233                pcldtaupi(i,k) = 0.0                 
234c
235                ELSE
236c                 
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
244c
245c-- cloud optical thickness :
246c [for liquid clouds, traditional formula,
247c for ice clouds, Ebert & Curry (1992)]
248c                 
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)
253c
254                ENDIF
255c
256          ENDDO           
257        ENDDO
258c
259      ELSE !--not ok_cdnc
260c
261c-prescribed cloud droplet radius
262c
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
275c
276c--computation of cloud optical depth and emissivity 
277c--in the general case
278c
279       DO k = 1, klev
280          DO i = 1, klon
281c
282             IF (pclc(i,k) .LE. seuil_neb) THEN
283c               
284c effective cloud droplet radius (microns) for liquid water clouds:
285c For output diagnostics cloud droplet effective radius [um]
286c we multiply here with f * xl (fraction of liquid water
287c clouds in the grid cell) to avoid problems in the averaging of the output.
288c In the output of IOIPSL, derive the REAL cloud droplet
289c effective radius as re/fl
290c
291                fl(i,k) = seuil_neb*(1.-zfice(i,k))           
292                re(i,k) = rad_chaud(i,k)*fl(i,k)
293                pclc(i,k)   = 0.0
294                pcltau(i,k) = 0.0
295                pclemi(i,k) = 0.0
296c
297             ELSE
298c
299c     -- liquid/ice cloud water paths:
300                 
301                zflwp_var= 1000.*(1.-zfice(i,k))*pqlwp(i,k)/pclc(i,k)
302     &               *rhodz(i,k)
303                zfiwp_var= 1000.*zfice(i,k)*pqlwp(i,k)/pclc(i,k)
304     &               *rhodz(i,k)
305c               
306c effective cloud droplet radius (microns) for liquid water clouds:
307c For output diagnostics cloud droplet effective radius [um]
308c we multiply here with f * xl (fraction of liquid water
309c clouds in the grid cell) to avoid problems in the averaging of the output.
310c In the output of IOIPSL, derive the REAL cloud droplet
311c effective radius as re/fl
312c
313                fl(i,k) = pclc(i,k)*(1.-zfice(i,k))           
314                re(i,k) = rad_chaud(i,k)*fl(i,k)
315c
316                rel = rad_chaud(i,k)
317c
318c for ice clouds: as a function of the ambiant temperature
319c [formula used by Iacobellis and Somerville (2000), with an
320c asymptotical value of 3.5 microns at T<-81.4 C added to be
321c consistent with observations of Heymsfield et al. 1986]:
322c 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as rei_max=61.29
323c
324                tc = t(i,k)-273.15
325                rei = d_rei_dt*tc + rei_max
326                if (tc.le.-81.4) rei = rei_min
327c
328c-- cloud optical thickness :
329c [for liquid clouds, traditional formula,
330c for ice clouds, Ebert & Curry (1992)]
331c                 
332                 if (zflwp_var.eq.0.) rel = 1.
333                 if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
334                 pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel )
335     &                 + zfiwp_var * (3.448e-03  + 2.431/rei)
336c
337c     -- cloud infrared emissivity:
338c     [the broadband infrared absorption coefficient is PARAMETERized
339c     as a function of the effective cld droplet radius]
340c     Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
341c
342                  k_ice = k_ice0 + 1.0/rei
343c                 
344                  pclemi(i,k) = 1.0
345     &                 - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var)
346c
347             ENDIF
348c
349             reliq(i,k)=rel
350             reice(i,k)=rei
351c
352             xflwp(i) = xflwp(i)+ xflwc(i,k) * rhodz(i,k)
353             xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * rhodz(i,k)
354c
355          ENDDO
356       ENDDO
357c
358c--if cloud droplet radius is fixed, then pcldtaupi=pcltau
359c
360      IF (.NOT.ok_cdnc) THEN
361         DO k = 1, klev
362            DO i = 1, klon
363               pcldtaupi(i,k)=pcltau(i,k)
364            ENDDO
365         ENDDO               
366      ENDIF
367C     
368C     COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
369c     IM cf. CR:test: calcul prenant ou non en compte le recouvrement
370c     initialisations
371c
372      DO i=1,klon
373         zclear(i)=1.
374         zcloud(i)=0.
375         zcloudh(i)=0.
376         zcloudm(i)=0.
377         zcloudl(i)=0.
378         pch(i)=1.0
379         pcm(i) = 1.0
380         pcl(i) = 1.0
381         pctlwp(i) = 0.0
382      ENDDO
383C
384c--calculation of liquid water path
385c
386      DO k = klev, 1, -1
387         DO i = 1, klon
388            pctlwp(i) = pctlwp(i)+ pqlwp(i,k)*rhodz(i,k)
389         ENDDO
390      ENDDO
391c     
392c--calculation of cloud properties with cloud overlap
393c
394      IF (NOVLP.EQ.1) THEN
395         DO k = klev, 1, -1
396            DO i = 1, klon
397               zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
398     &              /(1.-MIN(REAL(zcloud(i), kind=8),1.-ZEPSEC))
399               pct(i)=1.-zclear(i)
400               IF (paprs(i,k).LT.prmhc) THEN
401                  pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloudh(i)))
402     &                 /(1.-MIN(REAL(zcloudh(i), kind=8),1.-ZEPSEC))
403                  zcloudh(i)=pclc(i,k)
404               ELSE IF (paprs(i,k).GE.prmhc .AND.
405     &                 paprs(i,k).LT.prlmc) THEN
406                  pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloudm(i)))
407     &                 /(1.-MIN(REAL(zcloudm(i), kind=8),1.-ZEPSEC))
408                  zcloudm(i)=pclc(i,k)
409               ELSE IF (paprs(i,k).GE.prlmc) THEN
410                  pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloudl(i)))
411     &                 /(1.-MIN(REAL(zcloudl(i), kind=8),1.-ZEPSEC))
412                  zcloudl(i)=pclc(i,k)
413               endif
414               zcloud(i)=pclc(i,k)
415            ENDDO
416         ENDDO
417      ELSE IF (NOVLP.EQ.2) THEN
418         DO k = klev, 1, -1
419            DO i = 1, klon
420               zcloud(i)=MAX(pclc(i,k),zcloud(i))
421               pct(i)=zcloud(i)
422               IF (paprs(i,k).LT.prmhc) THEN
423                  pch(i) = MIN(pclc(i,k),pch(i))
424               ELSE IF (paprs(i,k).GE.prmhc .AND.
425     &                 paprs(i,k).LT.prlmc) THEN
426                  pcm(i) = MIN(pclc(i,k),pcm(i))
427               ELSE IF (paprs(i,k).GE.prlmc) THEN
428                  pcl(i) = MIN(pclc(i,k),pcl(i))
429               endif
430            ENDDO
431         ENDDO
432      ELSE IF (NOVLP.EQ.3) THEN
433         DO k = klev, 1, -1
434            DO i = 1, klon
435               zclear(i)=zclear(i)*(1.-pclc(i,k))
436               pct(i)=1-zclear(i)
437               IF (paprs(i,k).LT.prmhc) THEN
438                  pch(i) = pch(i)*(1.0-pclc(i,k))
439               ELSE IF (paprs(i,k).GE.prmhc .AND.
440     &                 paprs(i,k).LT.prlmc) THEN
441                  pcm(i) = pcm(i)*(1.0-pclc(i,k))
442               ELSE IF (paprs(i,k).GE.prlmc) THEN
443                  pcl(i) = pcl(i)*(1.0-pclc(i,k))
444               endif
445            ENDDO
446         ENDDO
447      ENDIF
448C     
449      DO i = 1, klon
450         pch(i)=1.-pch(i)
451         pcm(i)=1.-pcm(i)
452         pcl(i)=1.-pcl(i)
453      ENDDO
454c
455c ========================================================
456c DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
457c ========================================================
458c change by Nicolas Yan (LSCE)
459c Cloud Droplet Number Concentration (CDNC) : 3D variable
460c Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
461c Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
462c Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
463c Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
464c
465         IF (ok_cdnc) THEN
466c
467            DO k = 1, klev
468               DO i = 1, klon
469                  phase3d(i,k)=1-zfice(i,k)
470                  IF (pclc(i,k) .LE. seuil_neb) THEN
471                     lcc3d(i,k)=seuil_neb*phase3d(i,k)
472                  ELSE
473                     lcc3d(i,k)=pclc(i,k)*phase3d(i,k)
474                  ENDIF
475                  scdnc(i,k)=lcc3d(i,k)*cdnc(i,k) ! m-3
476               ENDDO
477            ENDDO
478c
479            DO i=1,klon
480               lcc(i)=0.
481               reffclwtop(i)=0.
482               cldncl(i)=0.
483               IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i) = 1.
484               IF(MAXIMUM) tcc(i) = 0.
485            ENDDO
486c
487            DO i=1,klon
488               DO k=klev-1,1,-1 !From TOA down
489c
490            ! Test, if the cloud optical depth exceeds the necessary
491            ! threshold:
492
493             IF (pcltau(i,k).GT.thres_tau
494     .            .AND. pclc(i,k).GT.thres_neb) THEN
495
496                  IF (MAXIMUM) THEN
497                    IF (FIRST) THEN
498                       write(*,*)'Hypothese de recouvrement: MAXIMUM'
499                       FIRST=.FALSE.
500                    ENDIF
501                    flag_max= -1.
502                    ftmp(i) = MAX(tcc(i),pclc(i,k))
503                  ENDIF
504
505                  IF (RANDOM) THEN
506                    IF (FIRST) THEN
507                       write(*,*)'Hypothese de recouvrement: RANDOM'
508                       FIRST=.FALSE.
509                    ENDIF
510                    flag_max= 1.
511                    ftmp(i) = tcc(i) * (1-pclc(i,k))
512                  ENDIF
513
514                  IF (MAXIMUM_RANDOM) THEN
515                    IF (FIRST) THEN
516                       write(*,*)'Hypothese de recouvrement: MAXIMUM_
517     .                         RANDOM'
518                       FIRST=.FALSE.
519                    ENDIF
520                    flag_max= 1.
521                    ftmp(i) = tcc(i) *
522     .              (1. - MAX(pclc(i,k),pclc(i,k+1))) /
523     .              (1. - MIN(pclc(i,k+1),1.-thres_neb))
524                  ENDIF
525c Effective radius of cloud droplet at top of cloud (m)
526                  reffclwtop(i) = reffclwtop(i) + rad_chaud(i,k) *
527     .           1.0E-06 * phase3d(i,k) * ( tcc(i) - ftmp(i))*flag_max
528c CDNC at top of cloud (m-3)
529                  cldncl(i) = cldncl(i) + cdnc(i,k) * phase3d(i,k) *
530     .                 (tcc(i) - ftmp(i))*flag_max
531c Liquid Cloud Content at top of cloud
532                  lcc(i) = lcc(i) + phase3d(i,k) * (tcc(i)-ftmp(i))*
533     .                    flag_max
534c Total Cloud Content at top of cloud
535                  tcc(i)=ftmp(i)
536             
537          ENDIF ! is there a visible, not-too-small cloud? 
538          ENDDO ! loop over k
539c
540          IF (RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i)
541c
542         ENDDO ! loop over i
543
544!! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC  REFFCLWS)
545            DO i = 1, klon
546               DO k = 1, klev
547! Weight to be used for outputs: eau_liquide*couverture nuageuse
548                  lcc3dcon(i,k) =rnebcon(i,k)*phase3d(i,k)*clwcon(i,k)  ! eau liquide convective
549                  lcc3dstra(i,k)=pclc(i,k)*pqlwp(i,k)*phase3d(i,k)
550                  lcc3dstra(i,k)=lcc3dstra(i,k)-lcc3dcon(i,k)           ! eau liquide stratiforme
551                  lcc3dstra(i,k)=MAX(lcc3dstra(i,k),0.0)
552! Compute cloud droplet radius as above in meter
553                  radius=1.1*((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
554     &               /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
555                  radius=MAX(radius, 5.e-6)
556! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
557                  reffclwc(i,k)=radius
558                  reffclwc(i,k)=reffclwc(i,k)*lcc3dcon(i,k)
559! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
560                  reffclws(i,k)=radius
561                  reffclws(i,k)=reffclws(i,k)*lcc3dstra(i,k)
562               ENDDO !klev
563            ENDDO !klon
564c
565c Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
566c
567            DO i = 1, klon
568               cldnvi(i)=0.
569               lcc_integrat(i)=0.
570               height(i)=0.
571               DO k = 1, klev
572                  cldnvi(i)=cldnvi(i)+cdnc(i,k)*lcc3d(i,k)*dh(i,k)
573                  lcc_integrat(i)=lcc_integrat(i)+lcc3d(i,k)*dh(i,k)
574                  height(i)=height(i)+dh(i,k)
575               ENDDO ! klev
576               lcc_integrat(i)=lcc_integrat(i)/height(i)
577               IF (lcc_integrat(i) .LE. 1.0E-03) THEN
578                  cldnvi(i)=cldnvi(i)*lcc(i)/seuil_neb
579               ELSE
580                  cldnvi(i)=cldnvi(i)*lcc(i)/lcc_integrat(i)
581               ENDIF
582            ENDDO ! klon
583           
584            DO i = 1, klon
585               DO k = 1, klev
586                IF (scdnc(i,k) .LE. 0.0)     scdnc(i,k)=0.0
587                IF (reffclws(i,k) .LE. 0.0)  reffclws(i,k)=0.0
588                IF (reffclwc(i,k) .LE. 0.0)  reffclwc(i,k)=0.0
589                IF (lcc3d(i,k) .LE. 0.0)     lcc3d(i,k)=0.0
590                IF (lcc3dcon(i,k) .LE. 0.0)  lcc3dcon(i,k)=0.0
591                IF (lcc3dstra(i,k) .LE. 0.0) lcc3dstra(i,k)=0.0
592               ENDDO
593               IF (reffclwtop(i) .LE. 0.0)  reffclwtop(i)=0.0
594               IF (cldncl(i) .LE. 0.0)      cldncl(i)=0.0
595               IF (cldnvi(i) .LE. 0.0)      cldnvi(i)=0.0
596               IF (lcc(i) .LE. 0.0)         lcc(i)=0.0
597            ENDDO
598c
599      ENDIF !ok_cdnc
600c
601      RETURN
602c
603      END
Note: See TracBrowser for help on using the repository browser.