source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/newmicro.F90 @ 3282

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

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