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

Last change on this file since 3138 was 3124, checked in by musat, 7 years ago

Add CMIP6' variables cldicemxrat, cldwatmxrat
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.9 KB
Line 
1! $Id: newmicro.F90 3124 2017-12-13 14:24:16Z jghattas $
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, icc3dcon, icc3dstra, zfice
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 rad_chaud(klon, klev) !--rayon pour les nuages chauds
155  REAL rad_chaud_pi(klon, klev) !--rayon pour les nuages chauds pre-industriels
156  REAL zflwp_var, zfiwp_var
157  REAL d_rei_dt
158
159  ! Abderrahmane oct 2009
160  REAL reliq(klon, klev), reice(klon, klev)
161  REAL reliq_pi(klon, klev), reice_pi(klon, klev)
162
163  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
164  ! FH : 2011/05/24
165
166  ! rei = ( rei_max - rei_min ) * T(°C) / 81.4 + rei_max
167  ! to be used for a temperature in celcius T(°C) < 0
168  ! rei=rei_min for T(°C) < -81.4
169
170  ! Calcul de la pente de la relation entre rayon effective des cristaux
171  ! et la température.
172  ! Pour retrouver les résultats numériques de la version d'origine,
173  ! on impose 0.71 quand on est proche de 0.71
174
175  d_rei_dt = (rei_max-rei_min)/81.4
176  IF (abs(d_rei_dt-0.71)<1.E-4) d_rei_dt = 0.71
177  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
178
179  ! Calculer l'epaisseur optique et l'emmissivite des nuages
180  ! IM inversion des DO
181
182  xflwp = 0.D0
183  xfiwp = 0.D0
184  xflwc = 0.D0
185  xfiwc = 0.D0
186
187  reliq = 0.
188  reice = 0.
189  reliq_pi = 0.
190  reice_pi = 0.
191
192  IF (iflag_t_glace.EQ.0) THEN
193    DO k = 1, klev
194      DO i = 1, klon
195        ! -layer calculation
196        rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
197        zrho(i, k) = pplay(i, k)/t(i, k)/rd ! kg/m3
198        dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
199        ! -Fraction of ice in cloud using a linear transition
200        zfice(i, k) = 1.0 - (t(i,k)-t_glace_min_old)/(t_glace_max_old-t_glace_min_old)
201        zfice(i, k) = min(max(zfice(i,k),0.0), 1.0)
202        ! -IM Total Liquid/Ice water content
203        xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k)
204        xfiwc(i, k) = zfice(i, k)*pqlwp(i, k)
205      END DO
206    END DO
207  ELSE ! of IF (iflag_t_glace.EQ.0)
208    DO k = 1, klev
209        CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:,k))
210 
211
212        ! JBM: icefrac_lsc is now contained icefrac_lsc_mod
213!       zfice(i, k) = icefrac_lsc(t(i,k), t_glace_min, &
214!                                 t_glace_max, exposant_glace)
215      DO i = 1, klon
216        ! -layer calculation
217        rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
218        zrho(i, k) = pplay(i, k)/t(i, k)/rd ! kg/m3
219        dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
220        ! -IM Total Liquid/Ice water content
221        xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k)
222        xfiwc(i, k) = zfice(i, k)*pqlwp(i, k)
223      END DO
224    END DO
225  ENDIF
226
227  IF (ok_cdnc) THEN
228
229    ! --we compute cloud properties as a function of the aerosol load
230
231    DO k = 1, klev
232      DO i = 1, klon
233
234        ! Formula "D" of Boucher and Lohmann, Tellus, 1995
235        ! Cloud droplet number concentration (CDNC) is restricted
236        ! to be within [20, 1000 cm^3]
237
238        ! --present-day case
239        cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
240          1.E-4))/log(10.))*1.E6 !-m-3
241        cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k)))
242
243        ! --pre-industrial case
244        cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
245          1.E-4))/log(10.))*1.E6 !-m-3
246        cdnc_pi(i, k) = min(1000.E6, max(20.E6,cdnc_pi(i,k)))
247
248        ! --present-day case
249        rad_chaud(i, k) = 1.1*((pqlwp(i,k)*pplay(i, &
250          k)/(rd*t(i,k)))/(4./3*rpi*1000.*cdnc(i,k)))**(1./3.)
251        rad_chaud(i, k) = max(rad_chaud(i,k)*1.E6, 5.)
252
253        ! --pre-industrial case
254        rad_chaud_pi(i, k) = 1.1*((pqlwp(i,k)*pplay(i, &
255          k)/(rd*t(i,k)))/(4./3.*rpi*1000.*cdnc_pi(i,k)))**(1./3.)
256        rad_chaud_pi(i, k) = max(rad_chaud_pi(i,k)*1.E6, 5.)
257
258        ! --pre-industrial case
259        ! --liquid/ice cloud water paths:
260        IF (pclc(i,k)<=seuil_neb) THEN
261
262          pcldtaupi(i, k) = 0.0
263
264        ELSE
265
266          zflwp_var = 1000.*(1.-zfice(i,k))*pqlwp(i, k)/pclc(i, k)* &
267            rhodz(i, k)
268          zfiwp_var = 1000.*zfice(i, k)*pqlwp(i, k)/pclc(i, k)*rhodz(i, k)
269          tc = t(i, k) - 273.15
270          rei = d_rei_dt*tc + rei_max
271          IF (tc<=-81.4) rei = rei_min
272
273          ! -- cloud optical thickness :
274          ! [for liquid clouds, traditional formula,
275          ! for ice clouds, Ebert & Curry (1992)]
276
277          IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
278          pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
279            zfiwp_var*(3.448E-03+2.431/rei)
280
281        END IF
282
283      END DO
284    END DO
285
286  ELSE !--not ok_cdnc
287
288    ! -prescribed cloud droplet radius
289
290    DO k = 1, min(3, klev)
291      DO i = 1, klon
292        rad_chaud(i, k) = rad_chau2
293        rad_chaud_pi(i, k) = rad_chau2
294      END DO
295    END DO
296    DO k = min(3, klev) + 1, klev
297      DO i = 1, klon
298        rad_chaud(i, k) = rad_chau1
299        rad_chaud_pi(i, k) = rad_chau1
300      END DO
301    END DO
302
303  END IF !--ok_cdnc
304
305  ! --computation of cloud optical depth and emissivity
306  ! --in the general case
307
308  DO k = 1, klev
309    DO i = 1, klon
310
311      IF (pclc(i,k)<=seuil_neb) THEN
312
313        ! effective cloud droplet radius (microns) for liquid water clouds:
314        ! For output diagnostics cloud droplet effective radius [um]
315        ! we multiply here with f * xl (fraction of liquid water
316        ! clouds in the grid cell) to avoid problems in the averaging of the
317        ! output.
318        ! In the output of IOIPSL, derive the REAL cloud droplet
319        ! effective radius as re/fl
320
321        fl(i, k) = seuil_neb*(1.-zfice(i,k))
322        re(i, k) = rad_chaud(i, k)*fl(i, k)
323        rel = 0.
324        rei = 0.
325        pclc(i, k) = 0.0
326        pcltau(i, k) = 0.0
327        pclemi(i, k) = 0.0
328
329      ELSE
330
331        ! -- liquid/ice cloud water paths:
332
333        zflwp_var = 1000.*(1.-zfice(i,k))*pqlwp(i, k)/pclc(i, k)*rhodz(i, k)
334        zfiwp_var = 1000.*zfice(i, k)*pqlwp(i, k)/pclc(i, k)*rhodz(i, k)
335
336        ! effective cloud droplet radius (microns) for liquid water clouds:
337        ! For output diagnostics cloud droplet effective radius [um]
338        ! we multiply here with f * xl (fraction of liquid water
339        ! clouds in the grid cell) to avoid problems in the averaging of the
340        ! output.
341        ! In the output of IOIPSL, derive the REAL cloud droplet
342        ! effective radius as re/fl
343
344        fl(i, k) = pclc(i, k)*(1.-zfice(i,k))
345        re(i, k) = rad_chaud(i, k)*fl(i, k)
346
347        rel = rad_chaud(i, k)
348
349        ! for ice clouds: as a function of the ambiant temperature
350        ! [formula used by Iacobellis and Somerville (2000), with an
351        ! asymptotical value of 3.5 microns at T<-81.4 C added to be
352        ! consistent with observations of Heymsfield et al. 1986]:
353        ! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
354        ! rei_max=61.29
355
356        tc = t(i, k) - 273.15
357        rei = d_rei_dt*tc + rei_max
358        IF (tc<=-81.4) rei = rei_min
359
360        ! -- cloud optical thickness :
361        ! [for liquid clouds, traditional formula,
362        ! for ice clouds, Ebert & Curry (1992)]
363
364        IF (zflwp_var==0.) rel = 1.
365        IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
366        pcltau(i, k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/ &
367          rei)
368
369        ! -- cloud infrared emissivity:
370        ! [the broadband infrared absorption coefficient is PARAMETERized
371        ! as a function of the effective cld droplet radius]
372        ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
373
374        k_ice = k_ice0 + 1.0/rei
375
376        pclemi(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
377
378      END IF
379
380      reice(i, k) = rei
381
382      xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
383      xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
384
385    END DO
386  END DO
387
388  ! --if cloud droplet radius is fixed, then pcldtaupi=pcltau
389
390  IF (.NOT. ok_cdnc) THEN
391    DO k = 1, klev
392      DO i = 1, klon
393        pcldtaupi(i, k) = pcltau(i, k)
394        reice_pi(i, k) = reice(i, k)
395      END DO
396    END DO
397  END IF
398
399  DO k = 1, klev
400    DO i = 1, klon
401      reliq(i, k) = rad_chaud(i, k)
402      reliq_pi(i, k) = rad_chaud_pi(i, k)
403      reice_pi(i, k) = reice(i, k)
404    END DO
405  END DO
406
407  ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
408  ! IM cf. CR:test: calcul prenant ou non en compte le recouvrement
409  ! initialisations
410
411  DO i = 1, klon
412    zclear(i) = 1.
413    zcloud(i) = 0.
414    zcloudh(i) = 0.
415    zcloudm(i) = 0.
416    zcloudl(i) = 0.
417    pch(i) = 1.0
418    pcm(i) = 1.0
419    pcl(i) = 1.0
420    pctlwp(i) = 0.0
421  END DO
422
423  ! --calculation of liquid water path
424
425  DO k = klev, 1, -1
426    DO i = 1, klon
427      pctlwp(i) = pctlwp(i) + pqlwp(i, k)*rhodz(i, k)
428    END DO
429  END DO
430
431  ! --calculation of cloud properties with cloud overlap
432
433  IF (novlp==1) THEN
434    DO k = klev, 1, -1
435      DO i = 1, klon
436        zclear(i) = zclear(i)*(1.-max(pclc(i,k),zcloud(i)))/(1.-min(real( &
437          zcloud(i),kind=8),1.-zepsec))
438        pct(i) = 1. - zclear(i)
439        IF (paprs(i,k)<prmhc) THEN
440          pch(i) = pch(i)*(1.-max(pclc(i,k),zcloudh(i)))/(1.-min(real(zcloudh &
441            (i),kind=8),1.-zepsec))
442          zcloudh(i) = pclc(i, k)
443        ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
444          pcm(i) = pcm(i)*(1.-max(pclc(i,k),zcloudm(i)))/(1.-min(real(zcloudm &
445            (i),kind=8),1.-zepsec))
446          zcloudm(i) = pclc(i, k)
447        ELSE IF (paprs(i,k)>=prlmc) THEN
448          pcl(i) = pcl(i)*(1.-max(pclc(i,k),zcloudl(i)))/(1.-min(real(zcloudl &
449            (i),kind=8),1.-zepsec))
450          zcloudl(i) = pclc(i, k)
451        END IF
452        zcloud(i) = pclc(i, k)
453      END DO
454    END DO
455  ELSE IF (novlp==2) THEN
456    DO k = klev, 1, -1
457      DO i = 1, klon
458        zcloud(i) = max(pclc(i,k), zcloud(i))
459        pct(i) = zcloud(i)
460        IF (paprs(i,k)<prmhc) THEN
461          pch(i) = min(pclc(i,k), pch(i))
462        ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
463          pcm(i) = min(pclc(i,k), pcm(i))
464        ELSE IF (paprs(i,k)>=prlmc) THEN
465          pcl(i) = min(pclc(i,k), pcl(i))
466        END IF
467      END DO
468    END DO
469  ELSE IF (novlp==3) THEN
470    DO k = klev, 1, -1
471      DO i = 1, klon
472        zclear(i) = zclear(i)*(1.-pclc(i,k))
473        pct(i) = 1 - zclear(i)
474        IF (paprs(i,k)<prmhc) THEN
475          pch(i) = pch(i)*(1.0-pclc(i,k))
476        ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
477          pcm(i) = pcm(i)*(1.0-pclc(i,k))
478        ELSE IF (paprs(i,k)>=prlmc) THEN
479          pcl(i) = pcl(i)*(1.0-pclc(i,k))
480        END IF
481      END DO
482    END DO
483  END IF
484
485  DO i = 1, klon
486    pch(i) = 1. - pch(i)
487    pcm(i) = 1. - pcm(i)
488    pcl(i) = 1. - pcl(i)
489  END DO
490
491  ! ========================================================
492  ! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
493  ! ========================================================
494  ! change by Nicolas Yan (LSCE)
495  ! Cloud Droplet Number Concentration (CDNC) : 3D variable
496  ! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
497  ! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
498  ! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
499  ! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
500
501  IF (ok_cdnc) THEN
502
503    DO k = 1, klev
504      DO i = 1, klon
505        phase3d(i, k) = 1 - zfice(i, k)
506        IF (pclc(i,k)<=seuil_neb) THEN
507          lcc3d(i, k) = seuil_neb*phase3d(i, k)
508        ELSE
509          lcc3d(i, k) = pclc(i, k)*phase3d(i, k)
510        END IF
511        scdnc(i, k) = lcc3d(i, k)*cdnc(i, k) ! m-3
512      END DO
513    END DO
514
515    DO i = 1, klon
516      lcc(i) = 0.
517      reffclwtop(i) = 0.
518      cldncl(i) = 0.
519      IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1.
520      IF (novlp.EQ.2) tcc(i) = 0.
521    END DO
522
523    DO i = 1, klon
524      DO k = klev - 1, 1, -1 !From TOA down
525
526          ! Test, if the cloud optical depth exceeds the necessary
527          ! threshold:
528
529        IF (pcltau(i,k)>thres_tau .AND. pclc(i,k)>thres_neb) THEN
530
531          IF (novlp.EQ.2) THEN
532            IF (first) THEN
533              WRITE (*, *) 'Hypothese de recouvrement: MAXIMUM'
534              first = .FALSE.
535            END IF
536            flag_max = -1.
537            ftmp(i) = max(tcc(i), pclc(i,k))
538          END IF
539
540          IF (novlp.EQ.3) THEN
541            IF (first) THEN
542              WRITE (*, *) 'Hypothese de recouvrement: RANDOM'
543              first = .FALSE.
544            END IF
545            flag_max = 1.
546            ftmp(i) = tcc(i)*(1-pclc(i,k))
547          END IF
548
549          IF (novlp.EQ.1) THEN
550            IF (first) THEN
551              WRITE (*, *) 'Hypothese de recouvrement: MAXIMUM_ &
552                &                                             &
553                &                                          RANDOM'
554              first = .FALSE.
555            END IF
556            flag_max = 1.
557            ftmp(i) = tcc(i)*(1.-max(pclc(i,k),pclc(i,k+1)))/(1.-min(pclc(i, &
558              k+1),1.-thres_neb))
559          END IF
560          ! Effective radius of cloud droplet at top of cloud (m)
561          reffclwtop(i) = reffclwtop(i) + rad_chaud(i, k)*1.0E-06*phase3d(i, &
562            k)*(tcc(i)-ftmp(i))*flag_max
563          ! CDNC at top of cloud (m-3)
564          cldncl(i) = cldncl(i) + cdnc(i, k)*phase3d(i, k)*(tcc(i)-ftmp(i))* &
565            flag_max
566          ! Liquid Cloud Content at top of cloud
567          lcc(i) = lcc(i) + phase3d(i, k)*(tcc(i)-ftmp(i))*flag_max
568          ! Total Cloud Content at top of cloud
569          tcc(i) = ftmp(i)
570
571        END IF ! is there a visible, not-too-small cloud?
572      END DO ! loop over k
573
574      IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1. - tcc(i)
575
576    END DO ! loop over i
577
578    ! ! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC
579    ! REFFCLWS)
580    DO i = 1, klon
581      DO k = 1, klev
582        ! Weight to be used for outputs: eau_liquide*couverture nuageuse
583        lcc3dcon(i, k) = rnebcon(i, k)*phase3d(i, k)*clwcon(i, k) ! eau liquide convective
584        lcc3dstra(i, k) = pclc(i, k)*pqlwp(i, k)*phase3d(i, k)
585        lcc3dstra(i, k) = lcc3dstra(i, k) - lcc3dcon(i, k) ! eau liquide stratiforme
586        lcc3dstra(i, k) = max(lcc3dstra(i,k), 0.0)
587        !FC pour la glace (CAUSES)
588        icc3dcon(i, k) = rnebcon(i, k)*(1-phase3d(i, k))*clwcon(i, k) !  glace convective
589        icc3dstra(i, k)= pclc(i, k)*pqlwp(i, k)*(1-phase3d(i, k))
590        icc3dstra(i, k) = icc3dstra(i, k) - icc3dcon(i, k) ! glace stratiforme
591        icc3dstra(i, k) = max( icc3dstra(i, k), 0.0)
592        !FC (CAUSES)
593
594        ! Compute cloud droplet radius as above in meter
595        radius = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3*rpi*1000.* &
596          cdnc(i,k)))**(1./3.)
597        radius = max(radius, 5.E-6)
598        ! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
599        reffclwc(i, k) = radius
600        reffclwc(i, k) = reffclwc(i, k)*lcc3dcon(i, k)
601        ! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
602        reffclws(i, k) = radius
603        reffclws(i, k) = reffclws(i, k)*lcc3dstra(i, k)
604      END DO !klev
605    END DO !klon
606
607    ! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
608
609    DO i = 1, klon
610      cldnvi(i) = 0.
611      lcc_integrat(i) = 0.
612      height(i) = 0.
613      DO k = 1, klev
614        cldnvi(i) = cldnvi(i) + cdnc(i, k)*lcc3d(i, k)*dh(i, k)
615        lcc_integrat(i) = lcc_integrat(i) + lcc3d(i, k)*dh(i, k)
616        height(i) = height(i) + dh(i, k)
617      END DO ! klev
618      lcc_integrat(i) = lcc_integrat(i)/height(i)
619      IF (lcc_integrat(i)<=1.0E-03) THEN
620        cldnvi(i) = cldnvi(i)*lcc(i)/seuil_neb
621      ELSE
622        cldnvi(i) = cldnvi(i)*lcc(i)/lcc_integrat(i)
623      END IF
624    END DO ! klon
625
626    DO i = 1, klon
627      DO k = 1, klev
628        IF (scdnc(i,k)<=0.0) scdnc(i, k) = 0.0
629        IF (reffclws(i,k)<=0.0) reffclws(i, k) = 0.0
630        IF (reffclwc(i,k)<=0.0) reffclwc(i, k) = 0.0
631        IF (lcc3d(i,k)<=0.0) lcc3d(i, k) = 0.0
632        IF (lcc3dcon(i,k)<=0.0) lcc3dcon(i, k) = 0.0
633        IF (lcc3dstra(i,k)<=0.0) lcc3dstra(i, k) = 0.0
634!FC (CAUSES)
635        IF (icc3dcon(i,k)<=0.0) icc3dcon(i, k) = 0.0
636        IF (icc3dstra(i,k)<=0.0) icc3dstra(i, k) = 0.0
637!FC (CAUSES)
638      END DO
639      IF (reffclwtop(i)<=0.0) reffclwtop(i) = 0.0
640      IF (cldncl(i)<=0.0) cldncl(i) = 0.0
641      IF (cldnvi(i)<=0.0) cldnvi(i) = 0.0
642      IF (lcc(i)<=0.0) lcc(i) = 0.0
643    END DO
644
645  END IF !ok_cdnc
646
647  RETURN
648
649END SUBROUTINE newmicro
Note: See TracBrowser for help on using the repository browser.