source: LMDZ6/trunk/libf/phylmd/newmicro.F90 @ 4553

Last change on this file since 4553 was 4535, checked in by evignon, 16 months ago

poursuite de la replay-isation de lscp en vue de la session
de reecriture de lscp_mod en juin

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