source: LMDZ6/trunk/libf/phylmd/lmdz_cloud_optics_prop.F90 @ 4687

Last change on this file since 4687 was 4683, checked in by evignon, 10 months ago

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