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

Last change on this file since 3999 was 3999, checked in by evignon, 3 years ago

commission de la nouvelle routine de condensation
grande echelle simplifiee (lscp, version epuree de fistilp)
et du schema de nuages de phase mixte (en developpement)

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