source: LMDZ6/trunk/libf/phylmd/lmdz_cloud_optics_prop.f90 @ 5496

Last change on this file since 5496 was 5496, checked in by evignon, 31 hours ago

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