source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phylmd/newmicro.F @ 5412

Last change on this file since 5412 was 1723, checked in by idelkadi, 12 years ago

Correction dans newmicro.F : initialisation des tableaux (rel et re)

  • 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 1723 2013-02-04 13:25:36Z evignon $
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                rel = 0.
294                rei = 0.
295                pclc(i,k)   = 0.0
296                pcltau(i,k) = 0.0
297                pclemi(i,k) = 0.0
298c
299             ELSE
300c
301c     -- 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)
307c               
308c effective cloud droplet radius (microns) for liquid water clouds:
309c For output diagnostics cloud droplet effective radius [um]
310c we multiply here with f * xl (fraction of liquid water
311c clouds in the grid cell) to avoid problems in the averaging of the output.
312c In the output of IOIPSL, derive the REAL cloud droplet
313c effective radius as re/fl
314c
315                fl(i,k) = pclc(i,k)*(1.-zfice(i,k))           
316                re(i,k) = rad_chaud(i,k)*fl(i,k)
317c
318                rel = rad_chaud(i,k)
319c
320c for ice clouds: as a function of the ambiant temperature
321c [formula used by Iacobellis and Somerville (2000), with an
322c asymptotical value of 3.5 microns at T<-81.4 C added to be
323c consistent with observations of Heymsfield et al. 1986]:
324c 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as rei_max=61.29
325c
326                tc = t(i,k)-273.15
327                rei = d_rei_dt*tc + rei_max
328                if (tc.le.-81.4) rei = rei_min
329c
330c-- cloud optical thickness :
331c [for liquid clouds, traditional formula,
332c for ice clouds, Ebert & Curry (1992)]
333c                 
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)
338c
339c     -- cloud infrared emissivity:
340c     [the broadband infrared absorption coefficient is PARAMETERized
341c     as a function of the effective cld droplet radius]
342c     Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
343c
344                  k_ice = k_ice0 + 1.0/rei
345c                 
346                  pclemi(i,k) = 1.0
347     &                 - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var)
348c
349             ENDIF
350c
351             reliq(i,k)=rel
352             reice(i,k)=rei
353c
354             xflwp(i) = xflwp(i)+ xflwc(i,k) * rhodz(i,k)
355             xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * rhodz(i,k)
356c
357          ENDDO
358       ENDDO
359c
360c--if cloud droplet radius is fixed, then pcldtaupi=pcltau
361c
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
369C     
370C     COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
371c     IM cf. CR:test: calcul prenant ou non en compte le recouvrement
372c     initialisations
373c
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
385C
386c--calculation of liquid water path
387c
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
393c     
394c--calculation of cloud properties with cloud overlap
395c
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
450C     
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
456c
457c ========================================================
458c DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
459c ========================================================
460c change by Nicolas Yan (LSCE)
461c Cloud Droplet Number Concentration (CDNC) : 3D variable
462c Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
463c Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
464c Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
465c Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
466c
467         IF (ok_cdnc) THEN
468c
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
480c
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
488c
489            DO i=1,klon
490               DO k=klev-1,1,-1 !From TOA down
491c
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(*,*)'Hypothese de recouvrement: MAXIMUM_
519     .                         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
527c 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
530c 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
533c Liquid Cloud Content at top of cloud
534                  lcc(i) = lcc(i) + phase3d(i,k) * (tcc(i)-ftmp(i))*
535     .                    flag_max
536c 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
541c
542          IF (RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i)
543c
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
566c
567c Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
568c
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
600c
601      ENDIF !ok_cdnc
602c
603      RETURN
604c
605      END
Note: See TracBrowser for help on using the repository browser.