source: LMDZ6/branches/Test_modipsl/libf/phylmd/newmicro.F90 @ 5456

Last change on this file since 5456 was 4562, checked in by evignon, 19 months ago

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