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

Last change on this file since 1779 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
RevLine 
[1279]1! $Id: newmicro.F 1723 2013-02-04 13:25:36Z fairhead $
[1523]2
3
[1279]4!     
[1712]5      SUBROUTINE newmicro (ok_cdnc, bl95_b0, bl95_b1,
6     .                  paprs, pplay,
[524]7     .                  t, pqlwp, pclc, pcltau, pclemi,
8     .                  pch, pcl, pcm, pct, pctlwp,
[1712]9     .                  xflwp, xfiwp, xflwc, xfiwc,
10     .                  mass_solu_aero, mass_solu_aero_pi,
11     .                  pcldtaupi, re, fl, reliq, reice)
12c
[766]13      USE dimphy
[1337]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
[524]18      IMPLICIT none
19c======================================================================
20c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
[1712]21c            O.   Boucher (LMD/CNRS) mise a jour en 201212
[524]22c Objet: Calculer epaisseur optique et emmissivite des nuages
23c======================================================================
24c Arguments:
[1712]25c ok_cdnc-input-L-flag pour calculer les rayons a partir des aerosols
26c
[524]27c t-------input-R-temperature
[1712]28c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere dans la partie nuageuse (kg/kg)
[524]29c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
[1279]30c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
[1712]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 (    -"-      )
[524]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           
[1712]40c
[524]41c pcltau--output-R-epaisseur optique des nuages
42c pclemi--output-R-emissivite des nuages (0 a 1)
[1712]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
[524]49c======================================================================
50C
51#include "YOMCST.h"
52#include "nuage.h"
[685]53#include "radepsi.h"
54#include "radopt.h"
[1712]55
[1337]56c choix de l'hypothese de recouvrememnt nuageuse
[1712]57      LOGICAL RANDOM, MAXIMUM_RANDOM, MAXIMUM
58      PARAMETER (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.)
59c
[1479]60      LOGICAL, SAVE :: FIRST=.TRUE.
61!$OMP THREADPRIVATE(FIRST)
[1337]62      INTEGER flag_max
[1712]63c
64c threshold PARAMETERs
[1337]65      REAL thres_tau,thres_neb
66      PARAMETER (thres_tau=0.3, thres_neb=0.001)
67c
[1712]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)
[524]73      REAL t(klon,klev)
74      REAL pclc(klon,klev)
75      REAL pqlwp(klon,klev)
[1712]76      REAL pcltau(klon,klev)
77      REAL pclemi(klon,klev)
78      REAL pcldtaupi(klon, klev)
[524]79c
[1712]80      REAL pct(klon)
81      REAL pcl(klon)
82      REAL pcm(klon)
83      REAL pch(klon)
84      REAL pctlwp(klon)
[524]85c
86      LOGICAL lo
87c
[1523]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
[1712]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
[1523]96      PARAMETER (prmhc = 440.*100., prlmc = 680.*100.)
[524]97C
98      INTEGER i, k
99      REAL xflwp(klon), xfiwp(klon)
100      REAL xflwc(klon,klev), xfiwc(klon,klev)
101c
[1712]102      REAL radius
103c
104      REAL coef_froi, coef_chau
[524]105      PARAMETER (coef_chau=0.13, coef_froi=0.09)
[1712]106c
[1286]107      REAL seuil_neb
108      PARAMETER (seuil_neb=0.001)
[1712]109c
[524]110      INTEGER nexpo ! exponentiel pour glace/eau
111      PARAMETER (nexpo=6)
[1712]112c      PARAMETER (nexpo=1)
[524]113
[1712]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
[524]119cjq for the aerosol indirect effect
120cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
121cjq     
[1279]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)
[524]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     
[1712]131      LOGICAL ok_cdnc
[524]132      REAL bl95_b0, bl95_b1     ! Parameter in B&L 95-Formula
133     
134cjq-end   
[685]135cIM cf. CR:parametres supplementaires
136      REAL zclear(klon)
137      REAL zcloud(klon)
[1523]138      REAL zcloudh(klon)
139      REAL zcloudm(klon)
140      REAL zcloudl(klon)
[1712]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
[1525]147      REAL d_rei_dt
[1146]148
[1279]149! Abderrahmane oct 2009
150      Real reliq(klon, klev), reice(klon, klev)
151
[1525]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
[1712]163c
[1525]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!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[524]167c
168c Calculer l'epaisseur optique et l'emmissivite des nuages
[1712]169c     IM inversion des DO
[524]170c
[1146]171      xflwp = 0.d0
172      xfiwp = 0.d0
173      xflwc = 0.d0
174      xfiwc = 0.d0
[1712]175c
[1334]176      reliq=0.
177      reice=0.
[1712]178c
[524]179      DO k = 1, klev
[1712]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)
[1146]192         ENDDO
193      ENDDO
[524]194
[1712]195      IF (ok_cdnc) THEN
196c
197c--we compute cloud properties as a function of the aerosol load
198c
199        DO k = 1, klev
[1146]200            DO i = 1, klon
[1712]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*
[1279]208     &               log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
[1712]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*
[1305]213     &               log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))
214     &               *1.e6 !-m-3
[1712]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
[524]273
[1712]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)
[1723]293                rel = 0.
294                rei = 0.
[1712]295                pclc(i,k)   = 0.0
296                pcltau(i,k) = 0.0
297                pclemi(i,k) = 0.0
298c
299             ELSE
300c
[1146]301c     -- liquid/ice cloud water paths:
302                 
[1712]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                 
[1279]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 )
[1146]337     &                 + zfiwp_var * (3.448e-03  + 2.431/rei)
[1712]338c
[1146]339c     -- cloud infrared emissivity:
[1712]340c     [the broadband infrared absorption coefficient is PARAMETERized
[1146]341c     as a function of the effective cld droplet radius]
342c     Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
[1712]343c
[1146]344                  k_ice = k_ice0 + 1.0/rei
[1712]345c                 
[1146]346                  pclemi(i,k) = 1.0
347     &                 - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var)
[1712]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
[1146]363         DO k = 1, klev
364            DO i = 1, klon
[1712]365               pcldtaupi(i,k)=pcltau(i,k)
[1146]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
[1712]373c
[685]374      DO i=1,klon
375         zclear(i)=1.
376         zcloud(i)=0.
[1523]377         zcloudh(i)=0.
378         zcloudm(i)=0.
379         zcloudl(i)=0.
[524]380         pch(i)=1.0
381         pcm(i) = 1.0
382         pcl(i) = 1.0
383         pctlwp(i) = 0.0
384      ENDDO
385C
[1712]386c--calculation of liquid water path
387c
[524]388      DO k = klev, 1, -1
[1146]389         DO i = 1, klon
[1712]390            pctlwp(i) = pctlwp(i)+ pqlwp(i,k)*rhodz(i,k)
[1146]391         ENDDO
392      ENDDO
[1712]393c     
394c--calculation of cloud properties with cloud overlap
395c
[1146]396      IF (NOVLP.EQ.1) THEN
397         DO k = klev, 1, -1
398            DO i = 1, klon
[685]399               zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
[1712]400     &              /(1.-MIN(REAL(zcloud(i), kind=8),1.-ZEPSEC))
[685]401               pct(i)=1.-zclear(i)
[1523]402               IF (paprs(i,k).LT.prmhc) THEN
403                  pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloudh(i)))
[1712]404     &                 /(1.-MIN(REAL(zcloudh(i), kind=8),1.-ZEPSEC))
[1523]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)))
[1712]409     &                 /(1.-MIN(REAL(zcloudm(i), kind=8),1.-ZEPSEC))
[1523]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)))
[1712]413     &                 /(1.-MIN(REAL(zcloudl(i), kind=8),1.-ZEPSEC))
[1523]414                  zcloudl(i)=pclc(i,k)
[685]415               endif
416               zcloud(i)=pclc(i,k)
[1146]417            ENDDO
418         ENDDO
419      ELSE IF (NOVLP.EQ.2) THEN
420         DO k = klev, 1, -1
421            DO i = 1, klon
[685]422               zcloud(i)=MAX(pclc(i,k),zcloud(i))
423               pct(i)=zcloud(i)
[1523]424               IF (paprs(i,k).LT.prmhc) THEN
[685]425                  pch(i) = MIN(pclc(i,k),pch(i))
[1523]426               ELSE IF (paprs(i,k).GE.prmhc .AND.
427     &                 paprs(i,k).LT.prlmc) THEN
[685]428                  pcm(i) = MIN(pclc(i,k),pcm(i))
[1523]429               ELSE IF (paprs(i,k).GE.prlmc) THEN
[685]430                  pcl(i) = MIN(pclc(i,k),pcl(i))
431               endif
[1146]432            ENDDO
433         ENDDO
434      ELSE IF (NOVLP.EQ.3) THEN
435         DO k = klev, 1, -1
436            DO i = 1, klon
[685]437               zclear(i)=zclear(i)*(1.-pclc(i,k))
438               pct(i)=1-zclear(i)
[1523]439               IF (paprs(i,k).LT.prmhc) THEN
[1146]440                  pch(i) = pch(i)*(1.0-pclc(i,k))
[1523]441               ELSE IF (paprs(i,k).GE.prmhc .AND.
442     &                 paprs(i,k).LT.prlmc) THEN
[1146]443                  pcm(i) = pcm(i)*(1.0-pclc(i,k))
[1523]444               ELSE IF (paprs(i,k).GE.prlmc) THEN
[1146]445                  pcl(i) = pcl(i)*(1.0-pclc(i,k))
[685]446               endif
[1146]447            ENDDO
[685]448         ENDDO
[1146]449      ENDIF
450C     
[524]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
[1712]456c
[1337]457c ========================================================
[1712]458c DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
[1337]459c ========================================================
[1712]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
[1337]469            DO k = 1, klev
470               DO i = 1, klon
[1712]471                  phase3d(i,k)=1-zfice(i,k)
[1337]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
[1712]480c
[1337]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
[1712]488c
[1337]489            DO i=1,klon
490               DO k=klev-1,1,-1 !From TOA down
[1712]491c
[1337]492            ! Test, if the cloud optical depth exceeds the necessary
493            ! threshold:
494
[1712]495             IF (pcltau(i,k).GT.thres_tau
496     .            .AND. pclc(i,k).GT.thres_neb) THEN
[1337]497
[1712]498                  IF (MAXIMUM) THEN
499                    IF (FIRST) THEN
[1337]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
[1712]507                  IF (RANDOM) THEN
508                    IF (FIRST) THEN
[1337]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
[1712]516                  IF (MAXIMUM_RANDOM) THEN
517                    IF (FIRST) THEN
[1337]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)
[1712]528                  reffclwtop(i) = reffclwtop(i) + rad_chaud(i,k) *
[1337]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
[1712]541c
542          IF (RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i)
543c
[1337]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
[1712]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)
[1337]558! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
[1712]559                  reffclwc(i,k)=radius
560                  reffclwc(i,k)=reffclwc(i,k)*lcc3dcon(i,k)
[1337]561! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
[1712]562                  reffclws(i,k)=radius
563                  reffclws(i,k)=reffclws(i,k)*lcc3dstra(i,k)
[1337]564               ENDDO !klev
565            ENDDO !klon
566c
[1712]567c Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
568c
[1337]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
[1712]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
[1337]593                IF (lcc3dstra(i,k) .LE. 0.0) lcc3dstra(i,k)=0.0
594               ENDDO
[1712]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
[1337]599            ENDDO
600c
[1712]601      ENDIF !ok_cdnc
602c
[524]603      RETURN
[1712]604c
[524]605      END
Note: See TracBrowser for help on using the repository browser.