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

Last change on this file since 4744 was 4715, checked in by Laurent Fairhead, 10 months ago

Final (hopefully) commit from the newmicro replayisation workshop. The final USE statements that were
still included in the cloud_optics_prop routine were moved to the call_cloud_optics_prop routine that
sets up the interface between LMDZ and the parametrization.
LF for LR, MCD, AI, EV, JBM
  • 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.1 KB
Line 
1! $Id: lmdz_cloud_optics_prop.F90 4715 2023-10-05 14:14:22Z jyg $
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.EQ.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 .EQ. 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 .EQ. 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 .GT. 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.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1.
603      IF (novlp.EQ.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.EQ.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.EQ.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.EQ.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.EQ.3 .OR. novlp.EQ.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  RETURN
733
734END SUBROUTINE cloud_optics_prop
735
736END MODULE lmdz_cloud_optics_prop
Note: See TracBrowser for help on using the repository browser.