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

Last change on this file since 4705 was 4704, checked in by evignon, 10 months ago

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