source: LMDZ6/branches/blowing_snow/libf/phylmd/newmicro.F90

Last change on this file was 4119, checked in by evignon, 3 years ago

add latitude_deg as input of newmicro in phylmdiso/physiq_mod

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