source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_cloud_optics_prop.F90 @ 5441

Last change on this file since 5441 was 5113, checked in by abarral, 5 months ago

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