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

Last change on this file since 3283 was 3283, checked in by musat, 6 years ago

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