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

Last change on this file since 4660 was 4639, checked in by evignon, 16 months ago

modification de Lea pour inclure une dependance de la phase a la temperature du sommet du nuage

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