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

Last change on this file since 4114 was 4114, checked in by evignon, 2 years ago

option pour inclure les effets radiatifs des chutes de neige
plus corrections du schema des nuages de phase mixte

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