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

Last change on this file since 5644 was 5643, checked in by aborella, 3 months ago

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