source: LMDZ6/trunk/libf/phylmd/newmicro.F90 @ 3102

Last change on this file since 3102 was 2596, checked in by musat, 8 years ago

In newmicro allow use of cloud overlap hypothesis defined
in radopt.h as for the radiation code.
In rrtm allow max and random cloud overlap hypothesis.
IM

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