source: LMDZ6/branches/Ocean_skin/libf/phylmd/newmicro.F90 @ 3811

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