source: LMDZ6/branches/LMDZ-QUEST/libf/phylmd/newmicro.F90 @ 5429

Last change on this file since 5429 was 3318, checked in by jghattas, 7 years ago

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