source: LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90 @ 5602

Last change on this file since 5602 was 5601, checked in by aborella, 2 months ago

Multiple changes:

  • tracers which were ratios are now absolute quantities. This is needed because when the ratio

is not defined, some aberrations may occur

  • added a new tracer for total specific humidity in contrails
  • rework of the mixing process for cirrus clouds (and contrails)
  • changed the numerical integration of ice crystals' sublimation
  • subroutines do not take real inputs anymore (at least klon tables)
  • added more radiative diagnostics for contrails
  • 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: 40.3 KB
Line 
1! $Id: lmdz_cloud_optics_prop.f90 5601 2025-04-02 14:04:40Z aborella $
2MODULE lmdz_cloud_optics_prop
3
4CONTAINS
5
6SUBROUTINE cloud_optics_prop(klon, klev, paprs, pplay, temp, radocond, picefra, pclc, &
7    pcltau, pclemi, pch, pcl, pcm, pct, radocondwp, xflwp, xfiwp, xflwc, xfiwc, &
8    mass_solu_aero, mass_solu_aero_pi, pcldtaupi, distcltop, temp_cltop, re, fl, reliq, reice, &
9    reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
10    reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
11    icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, &
12    !--AB contrails
13    contfra, qice_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
14    pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
15    xfiwp_nocont, xfiwc_nocont, reice_nocont)
16
17  USE lmdz_cloud_optics_prop_ini , ONLY : flag_aerosol, ok_cdnc
18  USE lmdz_cloud_optics_prop_ini , ONLY : lunout
19  USE lmdz_cloud_optics_prop_ini , ONLY : bl95_b0, bl95_b1
20  USE lmdz_cloud_optics_prop_ini , ONLY : latitude_deg
21  USE lmdz_cloud_optics_prop_ini , ONLY : iflag_t_glace
22  USE lmdz_cloud_optics_prop_ini , ONLY : cdnc_max, cdnc_max_m3
23  USE lmdz_cloud_optics_prop_ini , ONLY : cdnc_min, cdnc_min_m3
24  USE lmdz_cloud_optics_prop_ini , ONLY : thres_tau, thres_neb
25  USE lmdz_cloud_optics_prop_ini , ONLY : prmhc, prlmc
26  USE lmdz_cloud_optics_prop_ini , ONLY : coef_froi, coef_chau
27  USE lmdz_cloud_optics_prop_ini , ONLY : seuil_neb
28  USE lmdz_cloud_optics_prop_ini , ONLY : t_glace_min_old, t_glace_max_old
29  USE lmdz_cloud_optics_prop_ini , ONLY : k_ice0, df
30  USE lmdz_cloud_optics_prop_ini , ONLY : rg, rd, rpi
31  USE lmdz_cloud_optics_prop_ini , ONLY : rad_chau1, rad_chau2, rad_froid, iflag_rei
32  USE lmdz_cloud_optics_prop_ini , ONLY : ok_icefra_lscp, rei_max, rei_min
33  USE lmdz_cloud_optics_prop_ini , ONLY : rei_coef, rei_min_temp
34  USE lmdz_cloud_optics_prop_ini , ONLY : zepsec, novlp, iflag_ice_thermo, ok_new_lscp
35  USE lmdz_cloud_optics_prop_ini , ONLY : ok_plane_contrail, re_ice_crystals_contrails
36 
37
38
39
40  IMPLICIT NONE
41  ! ======================================================================
42  ! Authors: Z.X. Li (LMD/CNRS) date: 19930910
43  !          O.Boucher (LMD/CNRS) mise a jour en 201212
44  !          I. Musat (LMD/CNRS) : prise en compte de la meme hypothese
45  !                              de recouvrement pour les nuages que pour
46  !                              le rayonnement rrtm via le parametre
47  !                               novlp de radopt.h : 20160721
48  !          L.Fairheard, E.Vignon, JB Madeleine, L. Raillard, A. Idelkadi
49  !          M. Coulon-Decorzens: replayisation of the routine + cleaning
50  !                               and commentaries
51  !
52  ! Aim: compute condensate optical properties,
53  !      cloud optical depth and emissivity
54  ! ======================================================================
55 
56  ! List of arguments
57  !------------------
58
59  ! input:
60  INTEGER, INTENT(IN) :: klon, klev      ! number of horizontal and vertical grid points
61  REAL, INTENT(IN) :: paprs(klon, klev+1)! pressure at bottom interfaces [Pa]
62  REAL, INTENT(IN) :: pplay(klon, klev)  ! pressure at the middle of layers [Pa]
63  REAL, INTENT(IN) :: temp(klon, klev)   ! temperature [K]
64  REAL, INTENT(IN) :: radocond(klon, klev) ! cloud condensed water seen by radiation [kg/kg]
65  REAL, INTENT(IN) :: picefra(klon,klev) ! ice fraction in clouds from large scale condensation scheme [-]
66  REAL, INTENT(IN) :: rnebcon(klon,klev) ! convection cloud fraction [-]
67  REAL, INTENT(IN) :: ccwcon(klon,klev)  ! condensed water from deep convection [kg/kg]
68  ! jq for the aerosol indirect effect
69  ! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
70  REAL, INTENT(IN) :: mass_solu_aero(klon, klev)    ! total mass concentration for all soluble aerosols [ug m-3]
71  REAL, INTENT(IN) :: mass_solu_aero_pi(klon, klev) ! - (pre-industrial value)
72  REAL, INTENT(IN)  :: dNovrN(klon)         ! enhancement factor for cdnc
73  REAL, INTENT(OUT) :: distcltop(klon,klev) ! distance from large scale cloud top [m]
74  REAL, INTENT(OUT) :: temp_cltop(klon,klev)!temperature at large scale cloud top [K]
75
76  LOGICAL, INTENT(IN) :: ptconv(klon, klev) ! flag for grid points affected by deep convection
77
78  ! inout:
79  REAL, INTENT(INOUT) :: pclc(klon, klev) ! cloud fraction for radiation [-]
80
81  ! out:
82  REAL, INTENT(OUT) :: pct(klon)      ! 2D total cloud cover [-]
83  REAL, INTENT(OUT) :: pcl(klon)      ! 2D low cloud cover [-]
84  REAL, INTENT(OUT) :: pcm(klon)      ! 2D mid cloud cover [-]
85  REAL, INTENT(OUT) :: pch(klon)      ! 2D high cloud cover [-]
86  REAL, INTENT(OUT) :: radocondwp(klon) ! total condensed water path (seen by radiation) [kg/m2]
87  REAL, INTENT(OUT) :: xflwp(klon)    ! liquid water path (seen by radiation) [kg/m2]
88  REAL, INTENT(OUT) :: xfiwp(klon)    ! ice water path (seen by radiation) [kg/m2]
89  REAL, INTENT(OUT) :: xflwc(klon, klev) ! liquid water content seen by radiation [kg/kg]
90  REAL, INTENT(OUT) :: xfiwc(klon, klev) ! ice water content seen by radiation [kg/kg]
91  REAL, INTENT(OUT) :: re(klon, klev) ! cloud droplet effective radius multiplied by fl
92  REAL, INTENT(OUT) :: fl(klon, klev) ! xliq * rneb, denominator to re; fraction of liquid water clouds
93                                      ! introduced to avoid problems in the averaging of the output
94                                      ! water clouds within a grid cell
95
96  REAL, INTENT(OUT) :: pcltau(klon, klev) ! cloud optical depth [m]
97  REAL, INTENT(OUT) :: pclemi(klon, klev) ! cloud emissivity [-]
98  REAL, INTENT(OUT) :: pcldtaupi(klon, klev) ! pre-industrial value of cloud optical thickness, ie.
99                                             ! values of optical thickness that does not account
100                                             ! for aerosol effects on cloud droplet radius [m]
101
102  REAL, INTENT(OUT) :: reliq(klon, klev)   ! liquid droplet effective radius [m]
103  REAL, INTENT(OUT) :: reice(klon, klev)   ! ice effective radius [m]
104  REAL, INTENT(OUT) :: reliq_pi(klon, klev)! liquid droplet effective radius [m], pre-industrial
105  REAL, INTENT(OUT) :: reice_pi(klon, klev)! ice effective radius [m], pre-industrial
106  REAL, INTENT(OUT) :: scdnc(klon, klev)   ! cloud droplet number concentration, mean over the whole mesh [m-3]
107  REAL, INTENT(OUT) :: cldncl(klon)        ! cloud droplet number concentration at top of cloud [m-3]
108  REAL, INTENT(OUT) :: reffclwtop(klon)    ! effective radius of cloud droplet at top of cloud [m]
109  REAL, INTENT(OUT) :: lcc(klon)           ! liquid Cloud Content at top of cloud [kg/kg]
110  REAL, INTENT(OUT) :: reffclws(klon, klev)! stratiform cloud droplet effective radius
111  REAL, INTENT(OUT) :: reffclwc(klon, klev)! convective cloud droplet effective radius
112  REAL, INTENT(OUT) :: cldnvi(klon)        ! column Integrated cloud droplet Number [/m2]
113  REAL, INTENT(OUT) :: lcc3d(klon, klev)   ! cloud fraction for liquid part only [-]
114  REAL, INTENT(OUT) :: lcc3dcon(klon, klev)! cloud fraction for liquid part only, convective clouds [-]
115  REAL, INTENT(OUT) :: lcc3dstra(klon, klev)!cloud fraction for liquid part only, stratiform clouds [-]
116  REAL, INTENT(OUT) :: icc3dcon(klon, klev)! cloud fraction for liquid part only, convective clouds [-]
117  REAL, INTENT(OUT) :: icc3dstra(klon, klev)! cloud fraction for ice part only, stratiform clouds [-]
118  REAL, INTENT(INOUT) :: icefrac_optics(klon, klev)! ice fraction in clouds seen by radiation [-]
119
120  !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
121  REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
122  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! contrails condensed water [kg/kg]
123  REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
124  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
125  REAL, INTENT(OUT) :: xfiwp_nocont(klon)        ! ice water path (seen by radiation) without contrails [kg/m2]
126  REAL, INTENT(OUT) :: xfiwc_nocont(klon, klev)  ! ice water content seen by radiation without contrails [kg/kg]
127  REAL, INTENT(OUT) :: pclc_nocont(klon, klev)   ! cloud fraction for radiation without contrails [-]
128  REAL, INTENT(OUT) :: pcltau_nocont(klon, klev) ! cloud optical depth without contrails [-]
129  REAL, INTENT(OUT) :: pclemi_nocont(klon, klev) ! cloud emissivity without contrails [-]
130  REAL, INTENT(OUT) :: pcltau_cont(klon, klev)   ! contrails optical depth [-]
131  REAL, INTENT(OUT) :: pclemi_cont(klon, klev)   ! contrails emissivity [-]
132  REAL, INTENT(OUT) :: reice_nocont(klon, klev)  ! ice effective radius without contrails [microns]
133  !--AB
134
135  ! Local variables
136  !----------------
137
138  LOGICAL, SAVE :: first = .TRUE.
139  !$OMP THREADPRIVATE(first)
140  INTEGER flag_max
141
142  ! threshold PARAMETERs
143  REAL phase3d(klon, klev)
144  REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
145  LOGICAL lo
146  INTEGER i, k
147  REAL radius
148
149
150  REAL rel, tc, rei, iwc, dei, deimin, deimax
151  REAL k_ice
152
153  ! jq for the aerosol indirect effect
154  ! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
155  REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
156  REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
157  REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
158
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  REAL pclc_var
174
175
176  ! FH : 2011/05/24
177  ! rei = ( rei_max - rei_min ) * T(°C) / 81.4 + rei_max
178  ! to be used for a temperature in celcius T(°C) < 0
179  ! rei=rei_min for T(°C) < -81.4
180  ! Calcul de la pente de la relation entre rayon effective des cristaux
181  ! et la température Pour retrouver les résultats numériques de la version d'origine,
182  ! on impose 0.71 quand on est proche de 0.71
183  d_rei_dt = (rei_max-rei_min)/81.4
184  IF (abs(d_rei_dt-0.71)<1.E-4) d_rei_dt = 0.71
185
186  ! Calculer l'epaisseur optique et l'emmissivite des nuages
187  ! IM inversion des DO
188
189  xflwp = 0.D0
190  xfiwp = 0.D0
191  xflwc = 0.D0
192  xfiwc = 0.D0
193  !--AB
194  IF (ok_plane_contrail) THEN
195    xfiwp_nocont = 0.D0
196    xfiwc_nocont = 0.D0
197    reice_nocont = 0.
198  ENDIF
199
200  reliq = 0.
201  reice = 0.
202  reliq_pi = 0.
203  reice_pi = 0.
204
205  IF (iflag_t_glace.EQ.0) THEN
206    DO k = 1, klev
207      DO i = 1, klon
208        ! -layer calculation
209        rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
210        zrho(i, k) = pplay(i, k)/temp(i, k)/rd ! kg/m3
211        dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
212        ! -Fraction of ice in cloud using a linear transition
213        icefrac_optics(i, k) = 1.0 - (temp(i,k)-t_glace_min_old)/(t_glace_max_old-t_glace_min_old)
214        icefrac_optics(i, k) = min(max(icefrac_optics(i,k),0.0), 1.0)
215        ! -IM Total Liquid/Ice water content
216        xflwc(i, k) = (1.-icefrac_optics(i,k))*radocond(i, k)
217        xfiwc(i, k) = icefrac_optics(i, k)*radocond(i, k)
218      ENDDO
219    ENDDO
220  ELSE ! of IF (iflag_t_glace.EQ.0)
221    DO k = 1, klev
222
223
224!!$      IF (ok_new_lscp) THEN
225!!$          CALL icefrac_lscp(klon,temp(:,k),iflag_ice_thermo,distcltop(:,k),temp_cltop(:,k),icefrac_optics(:,k),dzfice(:,k))
226!!$      ELSE
227!!$          CALL icefrac_lsc(klon,temp(:,k),pplay(:,k)/paprs(:,1),icefrac_optics(:,k))
228!!$      ENDIF
229
230      DO i = 1, klon
231       
232        IF ((.NOT. ptconv(i,k)) .AND. ok_new_lscp .AND. ok_icefra_lscp) THEN
233        ! EV: take the ice fraction directly from the lscp code
234        ! consistent only for non convective grid points
235        ! critical for mixed phase clouds
236            icefrac_optics(i,k)=picefra(i,k)
237        ENDIF
238
239        ! -layer calculation
240        rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
241        zrho(i, k) = pplay(i, k)/temp(i, k)/rd ! kg/m3
242        dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
243        ! -IM Total Liquid/Ice water content
244        xflwc(i, k) = (1.-icefrac_optics(i,k))*radocond(i, k)
245        xfiwc(i, k) = icefrac_optics(i, k)*radocond(i, k)
246      ENDDO
247    ENDDO
248
249    !--AB
250    IF (ok_plane_contrail) THEN
251      !--If contrails are activated, we diagnose iwc without contrails
252      DO k = 1, klev
253        DO i = 1, klon
254          pclc_nocont(i,k) = pclc(i, k) - contfra(i, k)
255          xfiwc_nocont(i, k) = xfiwc(i, k) - qice_cont(i, k)
256        ENDDO
257      ENDDO
258    ENDIF
259  ENDIF
260
261
262
263
264
265
266  IF (ok_cdnc) THEN
267
268    ! --we compute cloud properties as a function of the aerosol load
269
270    DO k = 1, klev
271      DO i = 1, klon
272        ! Formula "D" of Boucher and Lohmann, Tellus, 1995
273        ! Cloud droplet number concentration (CDNC) is restricted
274        ! to be within [20, 1000 cm^3]
275
276        ! --pre-industrial case
277        cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
278          1.E-4))/log(10.))*1.E6 !-m-3
279        cdnc_pi(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc_pi(i,k)))
280
281      ENDDO
282    ENDDO
283
284    !--flag_aerosol=7 => MACv2SP climatology
285    !--in this case there is an enhancement factor
286    IF (flag_aerosol .EQ. 7) THEN
287
288      !--present-day
289      DO k = 1, klev
290        DO i = 1, klon
291          cdnc(i, k) = cdnc_pi(i,k)*dNovrN(i)
292        ENDDO
293      ENDDO
294
295    !--standard case
296    ELSE
297
298      DO k = 1, klev
299        DO i = 1, klon
300
301          ! Formula "D" of Boucher and Lohmann, Tellus, 1995
302          ! Cloud droplet number concentration (CDNC) is restricted
303          ! to be within [20, 1000 cm^3]
304
305          ! --present-day case
306          cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
307            1.E-4))/log(10.))*1.E6 !-m-3
308          cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k)))
309
310        ENDDO
311      ENDDO
312
313    ENDIF !--flag_aerosol
314
315    !--computing cloud droplet size
316    DO k = 1, klev
317      DO i = 1, klon
318
319        ! --present-day case
320        rad_chaud(i, k) = 1.1*((radocond(i,k)*pplay(i, &
321          k)/(rd*temp(i,k)))/(4./3*rpi*1000.*cdnc(i,k)))**(1./3.)
322        rad_chaud(i, k) = max(rad_chaud(i,k)*1.E6, 5.)
323
324        ! --pre-industrial case
325        rad_chaud_pi(i, k) = 1.1*((radocond(i,k)*pplay(i, &
326          k)/(rd*temp(i,k)))/(4./3.*rpi*1000.*cdnc_pi(i,k)))**(1./3.)
327        rad_chaud_pi(i, k) = max(rad_chaud_pi(i,k)*1.E6, 5.)
328
329        ! --pre-industrial case
330        ! --liquid/ice cloud water paths:
331        IF (pclc(i,k)<=seuil_neb) THEN
332
333          pcldtaupi(i, k) = 0.0
334
335        ELSE
336
337          zflwp_var = 1000.*(1.-icefrac_optics(i,k))*radocond(i, k)/pclc(i, k)* &
338            rhodz(i, k)
339          zfiwp_var = 1000.*icefrac_optics(i, k)*radocond(i, k)/pclc(i, k)*rhodz(i, k)
340          ! Calculation of ice cloud effective radius in micron
341          IF (iflag_rei .EQ. 2) THEN
342            ! in-cloud ice water content in g/m3
343            iwc = icefrac_optics(i,k) * radocond(i,k) / pclc(i,k) * zrho(i,k) * 1000.
344            ! this formula is a simplified version of the Sun 2001 one (as in the IFS model,
345            ! and which is activated for iflag_rei = 1).
346            ! In particular, the term in temperature**2 has been simplified.
347            ! The new coefs are tunable, and are by default set so that the results fit well
348            ! to the Sun 2001 formula
349            ! By default, rei_coef = 2.4 and rei_min_temp = 175.
350            ! The ranges of these parameters are determined so that the RMSE between this
351            ! formula and the one from Sun 2001 is less than 4 times the minimum RMSE
352            ! The ranges are [1.9, 2.9] for rei_coef and [160., 185.] for rei_min_temp
353            dei = rei_coef * (iwc**0.2445) * ( temp(i,k) - rei_min_temp )
354            ! we clip the results
355            !deimin = 20.
356            deimax = 155.
357            !dei = MIN(MAX(dei, deimin), deimax)
358            dei = MIN(dei, deimax)
359            ! formula to convert effective diameter in effective radius
360            rei = 3. * SQRT(3.) / 8. * dei
361            rei = MAX(rei, rei_min)
362          ELSEIF (iflag_rei .EQ. 1) THEN
363            ! when we account for precipitation in the radiation scheme,
364            ! It is recommended to use the rei formula from Sun and Rikkus 1999 with a revision
365            ! from Sun 2001 (as in the IFS model)
366            iwc=icefrac_optics(i, k)*radocond(i, k)/pclc(i,k)*zrho(i,k)*1000. !in cloud ice water content in g/m3
367            dei=(1.2351+0.0105*(temp(i,k)-273.15))*(45.8966*(iwc**0.2214) + &
368               & 0.7957*(iwc**0.2535)*(temp(i,k)-83.15))
369            !deimax=155.0
370            !deimin=20.+40*cos(abs(latitude_deg(i))/180.*RPI)
371            !Etienne: deimax and deimin controled by rei_max and rei_min in physiq.def
372            deimax=rei_max*2.0
373            deimin=2.0*rei_min+40*cos(abs(latitude_deg(i))/180.*RPI)
374            dei=min(dei,deimax)
375            dei=max(dei,deimin)
376            rei=3.*sqrt(3.)/8.*dei
377           ELSE
378            ! Default
379            ! for ice clouds: as a function of the ambiant temperature
380            ! [formula used by Iacobellis and Somerville (2000), with an
381            ! asymptotical value of 3.5 microns at T<-81.4 C added to be
382            ! consistent with observations of Heymsfield et al. 1986]:
383            ! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
384            ! rei_max=61.29
385            tc = temp(i, k) - 273.15
386            rei = d_rei_dt*tc + rei_max
387            IF (tc<=-81.4) rei = rei_min
388           ENDIF
389
390          ! -- cloud optical thickness :
391          ! [for liquid clouds, traditional formula,
392          ! for ice clouds, Ebert & Curry (1992)]
393
394          IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
395          IF ( ok_plane_contrail ) THEN
396            !--If contrails are activated, rei is a weighted average between the natural
397            !--rei and the contrails rei, with the weights being the fraction of natural
398            !--vs contrail cirrus in the gridbox
399            !--Beware, re_ice_crystals_contrails is in m, while rei is in microns
400            rei = ( rei * pclc_nocont(i,k) &
401                + re_ice_crystals_contrails * 1.E6 * contfra(i,k) ) / pclc(i,k)
402          ENDIF
403          pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
404            zfiwp_var*(3.448E-03+2.431/rei)
405
406        ENDIF
407
408      ENDDO
409    ENDDO
410
411  ELSE !--not ok_cdnc
412
413    ! -prescribed cloud droplet radius
414
415    DO k = 1, min(3, klev)
416      DO i = 1, klon
417        rad_chaud(i, k) = rad_chau2
418        rad_chaud_pi(i, k) = rad_chau2
419      ENDDO
420    ENDDO
421    DO k = min(3, klev) + 1, klev
422      DO i = 1, klon
423        rad_chaud(i, k) = rad_chau1
424        rad_chaud_pi(i, k) = rad_chau1
425      ENDDO
426    ENDDO
427
428  ENDIF !--ok_cdnc
429
430  ! --computation of cloud optical depth and emissivity
431  ! --in the general case
432
433  IF ( .NOT. ok_plane_contrail ) THEN
434
435  DO k = 1, klev
436    DO i = 1, klon
437
438      IF (pclc(i,k)<=seuil_neb) THEN
439
440        ! effective cloud droplet radius (microns) for liquid water clouds:
441        ! For output diagnostics cloud droplet effective radius [um]
442        ! we multiply here with f * xl (fraction of liquid water
443        ! clouds in the grid cell) to avoid problems in the averaging of the
444        ! output.
445        ! In the output of IOIPSL, derive the REAL cloud droplet
446        ! effective radius as re/fl
447
448        fl(i, k) = seuil_neb*(1.-icefrac_optics(i,k))
449        re(i, k) = rad_chaud(i, k)*fl(i, k)
450        rel = 0.
451        rei = 0.
452        pclc(i, k) = 0.0
453        pcltau(i, k) = 0.0
454        pclemi(i, k) = 0.0
455
456      ELSE
457
458        ! -- liquid/ice cloud water paths:
459
460        zflwp_var = 1000.*(1.-icefrac_optics(i,k))*radocond(i, k)/pclc(i, k)*rhodz(i, k)
461        zfiwp_var = 1000.*icefrac_optics(i, k)*radocond(i, k)/pclc(i, k)*rhodz(i, k)
462
463        ! effective cloud droplet radius (microns) for liquid water clouds:
464        ! For output diagnostics cloud droplet effective radius [um]
465        ! we multiply here with f Effective radius of cloud droplet at top of cloud (m)* xl (fraction of liquid water
466        ! clouds in the grid cell) to avoid problems in the averaging of the
467        ! output.
468        ! In the output of IOIPSL, derive the REAL cloud droplet
469        ! effective radius as re/fl
470
471        fl(i, k) = pclc(i, k)*(1.-icefrac_optics(i,k))
472        re(i, k) = rad_chaud(i, k)*fl(i, k)
473
474        rel = rad_chaud(i, k)
475
476        ! Calculation of ice cloud effective radius in micron
477
478        IF (iflag_rei .EQ. 2) THEN
479            ! in-cloud ice water content in g/m3
480            iwc = icefrac_optics(i,k) * radocond(i,k) / pclc(i,k) * zrho(i,k) * 1000.
481            ! this formula is a simplified version of the Sun 2001 one (as in the IFS model,
482            ! and which is activated for iflag_rei = 1).
483            ! In particular, the term in temperature**2 has been simplified.
484            ! The new coefs are tunable, and are by default set so that the results fit well
485            ! to the Sun 2001 formula
486            ! By default, rei_coef = 2.4 and rei_min_temp = 175.
487            ! The ranges of these parameters are determined so that the RMSE between this
488            ! formula and the one from Sun 2001 is less than 4 times the minimum RMSE
489            ! The ranges are [1.9, 2.9] for rei_coef and [160., 185.] for rei_min_temp
490            dei = rei_coef * (iwc**0.2445) * ( temp(i,k) - rei_min_temp )
491            ! we clip the results
492            !deimin = 20.
493            deimax = 155.
494            !dei = MIN(MAX(dei, deimin), deimax)
495            dei = MIN(dei, deimax)
496            ! formula to convert effective diameter to effective radius
497            rei = 3. * SQRT(3.) / 8. * dei
498            rei = MAX(rei, rei_min)
499           
500        ELSEIF (iflag_rei .EQ. 1) THEN
501            ! when we account for precipitation in the radiation scheme,
502            ! we use the rei formula from Sun and Rikkus 1999 with a revision
503            ! from Sun 2001 (as in the IFS model)
504            iwc=icefrac_optics(i, k)*radocond(i, k)/pclc(i,k)*zrho(i,k)*1000. !in cloud ice water content in g/m3
505            dei=(1.2351+0.0105*(temp(i,k)-273.15))*(45.8966*(iwc**0.2214) + &
506               &0.7957*(iwc**0.2535)*(temp(i,k)-83.15))
507            !deimax=155.0
508            !deimin=20.+40*cos(abs(latitude_deg(i))/180.*RPI)
509            !Etienne: deimax and deimin controled by rei_max and rei_min in physiq.def
510            deimax=rei_max*2.0
511            deimin=2.0*rei_min+40*cos(abs(latitude_deg(i))/180.*RPI)
512            dei=min(dei,deimax)
513            dei=max(dei,deimin)
514            rei=3.*sqrt(3.)/8.*dei
515       
516        ELSE
517            ! Default
518            ! for ice clouds: as a function of the ambiant temperature
519            ! [formula used by Iacobellis and Somerville (2000), with an
520            ! asymptotical value of 3.5 microns at T<-81.4 C added to be
521            ! consistent with observations of Heymsfield et al. 1986]:
522            ! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
523            ! rei_max=61.29
524            tc = temp(i, k) - 273.15
525            rei = d_rei_dt*tc + rei_max
526            IF (tc<=-81.4) rei = rei_min
527        ENDIF
528        ! -- cloud optical thickness :
529        ! [for liquid clouds, traditional formula,
530        ! for ice clouds, Ebert & Curry (1992)]
531
532        IF (zflwp_var==0.) rel = 1.
533        IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
534
535        pcltau(i, k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/ &
536          rei)
537
538        ! -- cloud infrared emissivity:
539        ! [the broadband infrared absorption coefficient is PARAMETERized
540        ! as a function of the effective cld droplet radius]
541        ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
542
543        k_ice = k_ice0 + 1.0/rei
544
545        pclemi(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
546
547      ENDIF
548
549      reice(i, k) = rei
550
551      xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
552      xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
553
554    ENDDO
555  ENDDO
556
557  ELSE !--AB if contrails are activated
558
559  DO k = 1, klev
560    DO i = 1, klon
561
562      IF (pclc(i,k)<=seuil_neb) THEN
563
564        ! effective cloud droplet radius (microns) for liquid water clouds:
565        ! For output diagnostics cloud droplet effective radius [um]
566        ! we multiply here with f * xl (fraction of liquid water
567        ! clouds in the grid cell) to avoid problems in the averaging of the
568        ! output.
569        ! In the output of IOIPSL, derive the REAL cloud droplet
570        ! effective radius as re/fl
571
572        fl(i, k) = seuil_neb*(1.-icefrac_optics(i,k))
573        re(i, k) = rad_chaud(i, k)*fl(i, k)
574        rel = 0.
575        rei = 0.
576        pclc(i, k) = 0.0
577        pcltau(i, k) = 0.0
578        pclemi(i, k) = 0.0
579
580        !--AB contrails
581        reice_nocont(i,k) = 0.
582        pclc_nocont(i,k) = 0.
583        pcltau_cont(i,k) = 0.
584        pclemi_cont(i,k) = 0.
585        pcltau_nocont(i,k) = 0.
586        pclemi_nocont(i,k) = 0.
587
588      ELSE
589
590        IF ( pclc_nocont(i,k) .GT. 0.01 * seuil_neb ) THEN
591          ! -- liquid/ice cloud water paths:
592
593          zflwp_var = 1000.*xflwc(i, k)/pclc(i, k)*rhodz(i, k)
594          zfiwp_var = 1000.*xfiwc_nocont(i, k)/pclc_nocont(i, k)*rhodz(i, k)
595
596          ! effective cloud droplet radius (microns) for liquid water clouds:
597          ! For output diagnostics cloud droplet effective radius [um]
598          ! we multiply here with f Effective radius of cloud droplet at top of cloud (m)* xl (fraction of liquid water
599          ! clouds in the grid cell) to avoid problems in the averaging of the
600          ! output.
601          ! In the output of IOIPSL, derive the REAL cloud droplet
602          ! effective radius as re/fl
603
604          fl(i, k) = pclc_nocont(i, k)*(1.-icefrac_optics(i,k))
605          re(i, k) = rad_chaud(i, k)*fl(i, k)
606
607          rel = rad_chaud(i, k)
608
609          ! Calculation of ice cloud effective radius in micron
610
611          IF (iflag_rei .EQ. 2) THEN
612              ! in-cloud ice water content in g/m3
613              iwc = xfiwc_nocont(i, k) / pclc_nocont(i,k) * zrho(i,k) * 1000.
614              ! this formula is a simplified version of the Sun 2001 one (as in the IFS model,
615              ! and which is activated for iflag_rei = 1).
616              ! In particular, the term in temperature**2 has been simplified.
617              ! The new coefs are tunable, and are by default set so that the results fit well
618              ! to the Sun 2001 formula
619              ! By default, rei_coef = 2.4 and rei_min_temp = 175.
620              ! The ranges of these parameters are determined so that the RMSE between this
621              ! formula and the one from Sun 2001 is less than 4 times the minimum RMSE
622              ! The ranges are [1.9, 2.9] for rei_coef and [160., 185.] for rei_min_temp
623              dei = rei_coef * (iwc**0.2445) * ( temp(i,k) - rei_min_temp )
624              ! we clip the results
625              !deimin = 20.
626              deimax = 155.
627              !dei = MIN(MAX(dei, deimin), deimax)
628              dei = MIN(dei, deimax)
629              ! formula to convert effective diameter to effective radius
630              rei = 3. * SQRT(3.) / 8. * dei
631              rei = MAX(rei, rei_min)
632             
633          ELSEIF (iflag_rei .EQ. 1) THEN
634              ! when we account for precipitation in the radiation scheme,
635              ! we use the rei formula from Sun and Rikkus 1999 with a revision
636              ! from Sun 2001 (as in the IFS model)
637              iwc = xfiwc_nocont(i, k) / pclc_nocont(i,k) * zrho(i,k) * 1000. !in cloud ice water content in g/m3
638              dei=(1.2351+0.0105*(temp(i,k)-273.15))*(45.8966*(iwc**0.2214) + &
639                 &0.7957*(iwc**0.2535)*(temp(i,k)-83.15))
640              !deimax=155.0
641              !deimin=20.+40*cos(abs(latitude_deg(i))/180.*RPI)
642              !Etienne: deimax and deimin controled by rei_max and rei_min in physiq.def
643              deimax=rei_max*2.0
644              deimin=2.0*rei_min+40*cos(abs(latitude_deg(i))/180.*RPI)
645              dei=min(dei,deimax)
646              dei=max(dei,deimin)
647              rei=3.*sqrt(3.)/8.*dei
648       
649          ELSE
650              ! Default
651              ! for ice clouds: as a function of the ambiant temperature
652              ! [formula used by Iacobellis and Somerville (2000), with an
653              ! asymptotical value of 3.5 microns at T<-81.4 C added to be
654              ! consistent with observations of Heymsfield et al. 1986]:
655              ! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
656              ! rei_max=61.29
657              tc = temp(i, k) - 273.15
658              rei = d_rei_dt*tc + rei_max
659              IF (tc<=-81.4) rei = rei_min
660          ENDIF
661          ! -- cloud optical thickness :
662          ! [for liquid clouds, traditional formula,
663          ! for ice clouds, Ebert & Curry (1992)]
664
665
666          IF (zflwp_var==0.) rel = 1.
667          IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
668
669          !--Diagnostics of clouds emissivity, optical depth and ice crystal radius
670          !--without contrails
671          reice_nocont(i,k) = rei
672          pcltau_nocont(i,k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/rei)
673          ! -- cloud infrared emissivity:
674          ! [the broadband infrared absorption coefficient is PARAMETERized
675          ! as a function of the effective cld droplet radius]
676          ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
677          k_ice = k_ice0 + 1.0/rei
678          pclemi_nocont(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
679
680        ELSE
681
682          rel = 1.
683          rei = 1.
684
685          !--Diagnostics of clouds emissivity, optical depth and ice crystal radius
686          !--without contrails
687          reice_nocont(i,k) = 1.
688          pclc_nocont(i,k) = 0.
689          pclc(i,k) = contfra(i,k)
690          pcltau_nocont(i, k) = 0.
691          pclemi_nocont(i, k) = 0.
692
693        ENDIF
694
695
696        IF ( contfra(i,k) .GT. 0.01 * seuil_neb ) THEN
697          !--Everything is the same but with contrails
698          zfiwp_var = 1000.*qice_cont(i, k)/contfra(i, k)*rhodz(i, k)
699          pclc_var = contfra(i,k)
700
701          pcltau_cont(i, k) = zfiwp_var*(3.448E-03+2.431/(re_ice_crystals_contrails*1.E6))
702          ! -- cloud infrared emissivity:
703          ! [the broadband infrared absorption coefficient is PARAMETERized
704          ! as a function of the effective cld droplet radius]
705          ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
706          k_ice = k_ice0 + 1.0/(re_ice_crystals_contrails*1.E6)
707          pclemi_cont(i, k) = 1.0 - exp(-df*k_ice*zfiwp_var)
708
709        ELSE
710          pclc_var = 0.
711          pcltau_cont(i, k) = 0.
712          pclemi_cont(i, k) = 0.
713        ENDIF
714
715        rei = ( re_ice_crystals_contrails*1.E6 * pclc_var &
716            + reice_nocont(i, k) * pclc_nocont(i,k) ) / ( pclc_var + pclc_nocont(i,k) )
717
718        zflwp_var = 1000.*xflwc(i, k)/pclc(i, k)*rhodz(i, k)
719        zfiwp_var = 1000.*xfiwc(i, k)/pclc(i, k)*rhodz(i, k)
720
721        pcltau(i,k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/rei)
722        ! -- cloud infrared emissivity:
723        ! [the broadband infrared absorption coefficient is PARAMETERized
724        ! as a function of the effective cld droplet radius]
725        ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
726        k_ice = k_ice0 + 1.0/rei
727        pclemi(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
728
729      ENDIF
730
731      reice(i, k) = rei
732
733      xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
734      xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
735      xfiwp_nocont(i) = xfiwp_nocont(i) + xfiwc_nocont(i, k)*rhodz(i, k)
736
737    ENDDO
738  ENDDO
739
740  ENDIF ! ok_plane_contrail
741
742  ! --if cloud droplet radius is fixed, then pcldtaupi=pcltau
743
744  IF (.NOT. ok_cdnc) THEN
745    DO k = 1, klev
746      DO i = 1, klon
747        pcldtaupi(i, k) = pcltau(i, k)
748        reice_pi(i, k) = reice(i, k)
749      ENDDO
750    ENDDO
751  ENDIF
752
753  DO k = 1, klev
754    DO i = 1, klon
755      reliq(i, k) = rad_chaud(i, k)
756      reliq_pi(i, k) = rad_chaud_pi(i, k)
757      reice_pi(i, k) = reice(i, k)
758    ENDDO
759  ENDDO
760
761  ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
762  ! IM cf. CR:test: calcul prenant ou non en compte le recouvrement
763  ! initialisations
764
765  DO i = 1, klon
766    zclear(i) = 1.
767    zcloud(i) = 0.
768    zcloudh(i) = 0.
769    zcloudm(i) = 0.
770    zcloudl(i) = 0.
771    pch(i) = 1.0
772    pcm(i) = 1.0
773    pcl(i) = 1.0
774    radocondwp(i) = 0.0
775  ENDDO
776
777  ! --calculation of liquid water path
778
779  DO k = klev, 1, -1
780    DO i = 1, klon
781      radocondwp(i) = radocondwp(i) + radocond(i, k)*rhodz(i, k)
782    ENDDO
783  ENDDO
784
785  ! --calculation of cloud properties with cloud overlap
786  ! choix de l'hypothese de recouvrement nuageuse via radopt.h (IM, 19.07.2016)
787  ! !novlp=1: max-random
788  ! !novlp=2: maximum
789  ! !novlp=3: random
790
791
792  IF (novlp==1) THEN
793    DO k = klev, 1, -1
794      DO i = 1, klon
795        zclear(i) = zclear(i)*(1.-max(pclc(i,k),zcloud(i)))/(1.-min(real( &
796          zcloud(i),kind=8),1.-zepsec))
797        pct(i) = 1. - zclear(i)
798        IF (paprs(i,k)<prmhc) THEN
799          pch(i) = pch(i)*(1.-max(pclc(i,k),zcloudh(i)))/(1.-min(real(zcloudh &
800            (i),kind=8),1.-zepsec))
801          zcloudh(i) = pclc(i, k)
802        ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
803          pcm(i) = pcm(i)*(1.-max(pclc(i,k),zcloudm(i)))/(1.-min(real(zcloudm &
804            (i),kind=8),1.-zepsec))
805          zcloudm(i) = pclc(i, k)
806        ELSE IF (paprs(i,k)>=prlmc) THEN
807          pcl(i) = pcl(i)*(1.-max(pclc(i,k),zcloudl(i)))/(1.-min(real(zcloudl &
808            (i),kind=8),1.-zepsec))
809          zcloudl(i) = pclc(i, k)
810        ENDIF
811        zcloud(i) = pclc(i, k)
812      ENDDO
813    ENDDO
814  ELSE IF (novlp==2) THEN
815    DO k = klev, 1, -1
816      DO i = 1, klon
817        zcloud(i) = max(pclc(i,k), zcloud(i))
818        pct(i) = zcloud(i)
819        IF (paprs(i,k)<prmhc) THEN
820          pch(i) = min(pclc(i,k), pch(i))
821        ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
822          pcm(i) = min(pclc(i,k), pcm(i))
823        ELSE IF (paprs(i,k)>=prlmc) THEN
824          pcl(i) = min(pclc(i,k), pcl(i))
825        ENDIF
826      ENDDO
827    ENDDO
828  ELSE IF (novlp==3) THEN
829    DO k = klev, 1, -1
830      DO i = 1, klon
831        zclear(i) = zclear(i)*(1.-pclc(i,k))
832        pct(i) = 1 - zclear(i)
833        IF (paprs(i,k)<prmhc) THEN
834          pch(i) = pch(i)*(1.0-pclc(i,k))
835        ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
836          pcm(i) = pcm(i)*(1.0-pclc(i,k))
837        ELSE IF (paprs(i,k)>=prlmc) THEN
838          pcl(i) = pcl(i)*(1.0-pclc(i,k))
839        ENDIF
840      ENDDO
841    ENDDO
842  ENDIF
843
844  DO i = 1, klon
845    pch(i) = 1. - pch(i)
846    pcm(i) = 1. - pcm(i)
847    pcl(i) = 1. - pcl(i)
848  ENDDO
849
850  !--AB we recompute high cloud cover with contrails only
851  !--and without contrails
852  IF ( ok_plane_contrail ) THEN
853    DO i = 1, klon
854      zclear(i) = 1.
855      zcloud(i) = 0.
856      zcloudh(i) = 0.
857      pch_nocont(i) = 1.
858      pct_cont(i) = 1.
859    ENDDO
860
861    DO k = klev, 1, -1
862      DO i = 1, klon
863        zclear(i) = zclear(i)*(1.-max(contfra(i,k),zcloud(i)))/(1.-min(real( &
864          zcloud(i),kind=8),1.-zepsec))
865        pct_cont(i) = 1. - zclear(i)
866        zcloud(i) = contfra(i,k)
867        IF (paprs(i,k)<prmhc) THEN
868          pch_nocont(i) = pch_nocont(i)*(1.-max(pclc_nocont(i,k),zcloudh(i)))/(1.-min(real( &
869            zcloudh(i),kind=8),1.-zepsec))
870          zcloudh(i) = pclc_nocont(i, k)
871        ENDIF
872      ENDDO
873    ENDDO
874
875    DO i = 1, klon
876      pch_nocont(i) = 1. - pch_nocont(i)
877    ENDDO
878  ENDIF ! ok_plane_contrail
879
880  ! ========================================================
881  ! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
882  ! ========================================================
883  ! change by Nicolas Yan (LSCE)
884  ! Cloud Droplet Number Concentration (CDNC) : 3D variable
885  ! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
886  ! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
887  ! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
888  ! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
889
890  IF (ok_cdnc) THEN
891
892    DO k = 1, klev
893      DO i = 1, klon
894        phase3d(i, k) = 1 - icefrac_optics(i, k)
895        IF (pclc(i,k)<=seuil_neb) THEN
896          lcc3d(i, k) = seuil_neb*phase3d(i, k)
897        ELSE
898          lcc3d(i, k) = pclc(i, k)*phase3d(i, k)
899        ENDIF
900        scdnc(i, k) = lcc3d(i, k)*cdnc(i, k) ! m-3
901      ENDDO
902    ENDDO
903
904    DO i = 1, klon
905      lcc(i) = 0.
906      reffclwtop(i) = 0.
907      cldncl(i) = 0.
908      IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1.
909      IF (novlp.EQ.2) tcc(i) = 0.
910    ENDDO
911
912    DO i = 1, klon
913      DO k = klev - 1, 1, -1 !From TOA down
914
915          ! Test, if the cloud optical depth exceeds the necessary
916          ! threshold:
917
918        IF (pcltau(i,k)>thres_tau .AND. pclc(i,k)>thres_neb) THEN
919
920          IF (novlp.EQ.2) THEN
921            IF (first) THEN
922              WRITE (*, *) 'Hypothese de recouvrement: MAXIMUM'
923              first = .FALSE.
924            ENDIF
925            flag_max = -1.
926            ftmp(i) = max(tcc(i), pclc(i,k))
927          ENDIF
928
929          IF (novlp.EQ.3) THEN
930            IF (first) THEN
931              WRITE (*, *) 'Hypothese de recouvrement: RANDOM'
932              first = .FALSE.
933            ENDIF
934            flag_max = 1.
935            ftmp(i) = tcc(i)*(1-pclc(i,k))
936          ENDIF
937
938          IF (novlp.EQ.1) THEN
939            IF (first) THEN
940              WRITE (*, *) 'Hypothese de recouvrement: MAXIMUM_ &
941                &                                             &
942                &                                          RANDOM'
943              first = .FALSE.
944            ENDIF
945            flag_max = 1.
946            ftmp(i) = tcc(i)*(1.-max(pclc(i,k),pclc(i,k+1)))/(1.-min(pclc(i, &
947              k+1),1.-thres_neb))
948          ENDIF
949          ! Effective radius of cloud droplet at top of cloud (m)
950          reffclwtop(i) = reffclwtop(i) + rad_chaud(i, k)*1.0E-06*phase3d(i, &
951            k)*(tcc(i)-ftmp(i))*flag_max
952          ! CDNC at top of cloud (m-3)
953          cldncl(i) = cldncl(i) + cdnc(i, k)*phase3d(i, k)*(tcc(i)-ftmp(i))* &
954            flag_max
955          ! Liquid Cloud Content at top of cloud
956          lcc(i) = lcc(i) + phase3d(i, k)*(tcc(i)-ftmp(i))*flag_max
957          ! Total Cloud Content at top of cloud
958          tcc(i) = ftmp(i)
959
960        ENDIF ! is there a visible, not-too-small cloud?
961      ENDDO ! loop over k
962
963      IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1. - tcc(i)
964
965    ENDDO ! loop over i
966
967    ! ! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC
968    ! REFFCLWS)
969    DO i = 1, klon
970      DO k = 1, klev
971        ! Weight to be used for outputs: eau_liquide*couverture nuageuse
972        lcc3dcon(i, k) = rnebcon(i, k)*phase3d(i, k)*ccwcon(i, k) ! eau liquide convective
973        lcc3dstra(i, k) = pclc(i, k)*radocond(i, k)*phase3d(i, k)
974        lcc3dstra(i, k) = lcc3dstra(i, k) - lcc3dcon(i, k) ! eau liquide stratiforme
975        lcc3dstra(i, k) = max(lcc3dstra(i,k), 0.0)
976        !FC pour la glace (CAUSES)
977        icc3dcon(i, k) = rnebcon(i, k)*(1-phase3d(i, k))*ccwcon(i, k) !  glace convective
978        icc3dstra(i, k)= pclc(i, k)*radocond(i, k)*(1-phase3d(i, k))
979        icc3dstra(i, k) = icc3dstra(i, k) - icc3dcon(i, k) ! glace stratiforme
980        icc3dstra(i, k) = max( icc3dstra(i, k), 0.0)
981        !FC (CAUSES)
982
983        ! Compute cloud droplet radius as above in meter
984        radius = 1.1*((radocond(i,k)*pplay(i,k)/(rd*temp(i,k)))/(4./3*rpi*1000.* &
985          cdnc(i,k)))**(1./3.)
986        radius = max(radius, 5.E-6)
987        ! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
988        reffclwc(i, k) = radius
989        reffclwc(i, k) = reffclwc(i, k)*lcc3dcon(i, k)
990        ! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
991        reffclws(i, k) = radius
992        reffclws(i, k) = reffclws(i, k)*lcc3dstra(i, k)
993      ENDDO !klev
994    ENDDO !klon
995
996    ! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
997
998    DO i = 1, klon
999      cldnvi(i) = 0.
1000      lcc_integrat(i) = 0.
1001      height(i) = 0.
1002      DO k = 1, klev
1003        cldnvi(i) = cldnvi(i) + cdnc(i, k)*lcc3d(i, k)*dh(i, k)
1004        lcc_integrat(i) = lcc_integrat(i) + lcc3d(i, k)*dh(i, k)
1005        height(i) = height(i) + dh(i, k)
1006      ENDDO ! klev
1007      lcc_integrat(i) = lcc_integrat(i)/height(i)
1008      IF (lcc_integrat(i)<=1.0E-03) THEN
1009        cldnvi(i) = cldnvi(i)*lcc(i)/seuil_neb
1010      ELSE
1011        cldnvi(i) = cldnvi(i)*lcc(i)/lcc_integrat(i)
1012      ENDIF
1013    ENDDO ! klon
1014
1015    DO i = 1, klon
1016      DO k = 1, klev
1017        IF (scdnc(i,k)<=0.0) scdnc(i, k) = 0.0
1018        IF (reffclws(i,k)<=0.0) reffclws(i, k) = 0.0
1019        IF (reffclwc(i,k)<=0.0) reffclwc(i, k) = 0.0
1020        IF (lcc3d(i,k)<=0.0) lcc3d(i, k) = 0.0
1021        IF (lcc3dcon(i,k)<=0.0) lcc3dcon(i, k) = 0.0
1022        IF (lcc3dstra(i,k)<=0.0) lcc3dstra(i, k) = 0.0
1023!FC (CAUSES)
1024        IF (icc3dcon(i,k)<=0.0) icc3dcon(i, k) = 0.0
1025        IF (icc3dstra(i,k)<=0.0) icc3dstra(i, k) = 0.0
1026!FC (CAUSES)
1027      ENDDO
1028      IF (reffclwtop(i)<=0.0) reffclwtop(i) = 0.0
1029      IF (cldncl(i)<=0.0) cldncl(i) = 0.0
1030      IF (cldnvi(i)<=0.0) cldnvi(i) = 0.0
1031      IF (lcc(i)<=0.0) lcc(i) = 0.0
1032    ENDDO
1033
1034  ENDIF !ok_cdnc
1035
1036  first=.false. !to be sure
1037
1038  RETURN
1039
1040END SUBROUTINE cloud_optics_prop
1041
1042END MODULE lmdz_cloud_optics_prop
Note: See TracBrowser for help on using the repository browser.