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

Last change on this file since 3274 was 3274, checked in by oboucher, 6 years ago

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