GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/newmicro.F90 Lines: 125 244 51.2 %
Date: 2023-06-30 12:51:15 Branches: 94 200 47.0 %

Line Branch Exec Source
1
! $Id: newmicro.F90 4562 2023-06-07 13:33:57Z evignon $
2
3
12159804
SUBROUTINE newmicro(flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, picefra, pclc, &
4
72
    pcltau, pclemi, pch, pcl, pcm, pct, pctlwp, xflwp, xfiwp, xflwc, xfiwc, &
5
    mass_solu_aero, mass_solu_aero_pi, pcldtaupi, latitude_deg,distcltop, 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
  ! distcltop ---input- distance from cloud top
45
46
  ! re------output-R-Cloud droplet effective radius multiplied by fl [um]
47
  ! fl------output-R-Denominator to re, introduced to avoid problems in
48
  ! the averaging of the output. fl is the fraction of liquid
49
  ! water clouds within a grid cell
50
51
  ! pcltau--output-R-epaisseur optique des nuages
52
  ! pclemi--output-R-emissivite des nuages (0 a 1)
53
  ! pcldtaupi-output-R-pre-industrial value of cloud optical thickness,
54
55
  ! pcl-output-R-2D low-level cloud cover
56
  ! pcm-output-R-2D mid-level cloud cover
57
  ! pch-output-R-2D high-level cloud cover
58
  ! pct-output-R-2D total cloud cover
59
  ! ======================================================================
60
61
  include "YOMCST.h"
62
  include "nuage.h"
63
  include "radepsi.h"
64
  include "radopt.h"
65
  include "clesphys.h"
66
67
  ! choix de l'hypothese de recouvrement nuageuse via radopt.h (IM, 19.07.2016)
68
  ! !novlp=1: max-random
69
  ! !novlp=2: maximum
70
  ! !novlp=3: random
71
! LOGICAL random, maximum_random, maximum
72
! PARAMETER (random=.FALSE., maximum_random=.TRUE., maximum=.FALSE.)
73
74
  LOGICAL, SAVE :: first = .TRUE.
75
  !$OMP THREADPRIVATE(FIRST)
76
  INTEGER flag_max
77
78
  ! threshold PARAMETERs
79
  REAL thres_tau, thres_neb
80
  PARAMETER (thres_tau=0.3, thres_neb=0.001)
81
82
144
  REAL phase3d(klon, klev)
83
144
  REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
84
85
  REAL paprs(klon, klev+1)
86
  REAL pplay(klon, klev)
87
  REAL t(klon, klev)
88
  REAL pclc(klon, klev)
89
  REAL pqlwp(klon, klev), picefra(klon,klev)
90
  REAL pcltau(klon, klev)
91
  REAL pclemi(klon, klev)
92
  REAL pcldtaupi(klon, klev)
93
  REAL latitude_deg(klon)
94
95
  REAL pct(klon)
96
  REAL pcl(klon)
97
  REAL pcm(klon)
98
  REAL pch(klon)
99
  REAL pctlwp(klon)
100
101
  REAL distcltop(klon,klev)
102
  LOGICAL lo
103
104
  ! !Abderr modif JL mail du 19.01.2011 18:31
105
  ! REAL cetahb, cetamb
106
  ! PARAMETER (cetahb = 0.45, cetamb = 0.80)
107
  ! Remplacer
108
  ! cetahb*paprs(i,1) par  prmhc
109
  ! cetamb*paprs(i,1) par  prlmc
110
  REAL prmhc ! Pressure between medium and high level cloud in Pa
111
  REAL prlmc ! Pressure between low and medium level cloud in Pa
112
  PARAMETER (prmhc=440.*100., prlmc=680.*100.)
113
114
  INTEGER i, k
115
  REAL xflwp(klon), xfiwp(klon)
116
  REAL xflwc(klon, klev), xfiwc(klon, klev)
117
118
  REAL radius
119
120
  REAL coef_froi, coef_chau
121
  PARAMETER (coef_chau=0.13, coef_froi=0.09)
122
123
  REAL seuil_neb
124
  PARAMETER (seuil_neb=0.001)
125
126
! JBM (3/14) nexpo is replaced by exposant_glace
127
! INTEGER nexpo ! exponentiel pour glace/eau
128
! PARAMETER (nexpo=6)
129
! PARAMETER (nexpo=1)
130
! if iflag_t_glace=0, the old values are used:
131
  REAL, PARAMETER :: t_glace_min_old = 258.
132
  REAL, PARAMETER :: t_glace_max_old = 273.13
133
134
  REAL rel, tc, rei, iwc, dei, deimin, deimax
135
  REAL k_ice0, k_ice, df
136
  PARAMETER (k_ice0=0.005) ! units=m2/g
137
  PARAMETER (df=1.66) ! diffusivity factor
138
139
  ! jq for the aerosol indirect effect
140
  ! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
141
  ! jq
142
  REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols [ug m-3]
143
  REAL mass_solu_aero_pi(klon, klev) ! - " - (pre-industrial value)
144
144
  REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
145
  REAL re(klon, klev) ! cloud droplet effective radius [um]
146
144
  REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
147
  REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
148
149
  REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds
150
  ! within the grid cell)
151
152
  INTEGER flag_aerosol
153
  LOGICAL ok_cdnc
154
  REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
155
156
  ! jq-end
157
  ! IM cf. CR:parametres supplementaires
158
144
  REAL dzfice(klon,klev)
159
144
  REAL zclear(klon)
160
144
  REAL zcloud(klon)
161
144
  REAL zcloudh(klon)
162
144
  REAL zcloudm(klon)
163
144
  REAL zcloudl(klon)
164
144
  REAL rhodz(klon, klev) !--rho*dz pour la couche
165
144
  REAL zrho(klon, klev) !--rho pour la couche
166
144
  REAL dh(klon, klev) !--dz pour la couche
167
144
  REAL rad_chaud(klon, klev) !--rayon pour les nuages chauds
168
144
  REAL rad_chaud_pi(klon, klev) !--rayon pour les nuages chauds pre-industriels
169
  REAL zflwp_var, zfiwp_var
170
  REAL d_rei_dt
171
172
  ! Abderrahmane oct 2009
173
  REAL reliq(klon, klev), reice(klon, klev)
174
  REAL reliq_pi(klon, klev), reice_pi(klon, klev)
175
176
  REAL,SAVE :: cdnc_min=-1.
177
  REAL,SAVE :: cdnc_min_m3
178
  !$OMP THREADPRIVATE(cdnc_min,cdnc_min_m3)
179
  REAL,SAVE :: cdnc_max=-1.
180
  REAL,SAVE :: cdnc_max_m3
181
  !$OMP THREADPRIVATE(cdnc_max,cdnc_max_m3)
182
183
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
184
  ! FH : 2011/05/24
185
186
  ! rei = ( rei_max - rei_min ) * T(°C) / 81.4 + rei_max
187
  ! to be used for a temperature in celcius T(°C) < 0
188
  ! rei=rei_min for T(°C) < -81.4
189
190
  ! Calcul de la pente de la relation entre rayon effective des cristaux
191
  ! et la température.
192
  ! Pour retrouver les résultats numériques de la version d'origine,
193
  ! on impose 0.71 quand on est proche de 0.71
194
195
72
  if (first) THEN
196
1
      call getin_p('cdnc_min',cdnc_min)
197
1
      cdnc_min_m3=cdnc_min*1.E6
198
1
      IF (cdnc_min_m3<0.) cdnc_min_m3=20.E6 ! astuce pour retrocompatibilite
199
1
      write(lunout,*)'cdnc_min=', cdnc_min_m3/1.E6
200
1
      call getin_p('cdnc_max',cdnc_max)
201
1
      cdnc_max_m3=cdnc_max*1.E6
202
1
      IF (cdnc_max_m3<0.) cdnc_max_m3=1000.E6 ! astuce pour retrocompatibilite
203
1
      write(lunout,*)'cdnc_max=', cdnc_max_m3/1.E6
204
  ENDIF
205
206
72
  d_rei_dt = (rei_max-rei_min)/81.4
207
72
  IF (abs(d_rei_dt-0.71)<1.E-4) d_rei_dt = 0.71
208
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
209
210
  ! Calculer l'epaisseur optique et l'emmissivite des nuages
211
  ! IM inversion des DO
212
213
71640
  xflwp = 0.D0
214
71640
  xfiwp = 0.D0
215

2794032
  xflwc = 0.D0
216

2794032
  xfiwc = 0.D0
217
218

2794032
  reliq = 0.
219

2794032
  reice = 0.
220

2794032
  reliq_pi = 0.
221

2794032
  reice_pi = 0.
222
223
72
  IF (iflag_t_glace.EQ.0) THEN
224
    DO k = 1, klev
225
      DO i = 1, klon
226
        ! -layer calculation
227
        rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
228
        zrho(i, k) = pplay(i, k)/t(i, k)/rd ! kg/m3
229
        dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
230
        ! -Fraction of ice in cloud using a linear transition
231
        zfice(i, k) = 1.0 - (t(i,k)-t_glace_min_old)/(t_glace_max_old-t_glace_min_old)
232
        zfice(i, k) = min(max(zfice(i,k),0.0), 1.0)
233
        ! -IM Total Liquid/Ice water content
234
        xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k)
235
        xfiwc(i, k) = zfice(i, k)*pqlwp(i, k)
236
      ENDDO
237
    ENDDO
238
  ELSE ! of IF (iflag_t_glace.EQ.0)
239
2880
    DO k = 1, klev
240
241
! JBM: icefrac_lsc is now contained icefrac_lsc_mod
242
!       zfice(i, k) = icefrac_lsc(t(i,k), t_glace_min, &
243
!                                 t_glace_max, exposant_glace)
244
245
2808
      IF (ok_new_lscp) THEN
246
2808
          CALL icefrac_lscp(klon,t(:,k),iflag_ice_thermo,distcltop(:,k),zfice(:,k),dzfice(:,k))
247
      ELSE
248
          CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:,k))
249
      ENDIF
250
251
2794032
      DO i = 1, klon
252
253

2791152
        IF ((.NOT. ptconv(i,k)) .AND. ok_new_lscp .AND. ok_icefra_lscp) THEN
254
        ! EV: take the ice fraction directly from the lscp code
255
        ! consistent only for non convective grid points
256
        ! critical for mixed phase clouds
257
2666973
            zfice(i,k)=picefra(i,k)
258
        ENDIF
259
260
        ! -layer calculation
261
2791152
        rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
262
2791152
        zrho(i, k) = pplay(i, k)/t(i, k)/rd ! kg/m3
263
2791152
        dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
264
        ! -IM Total Liquid/Ice water content
265
2791152
        xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k)
266
2793960
        xfiwc(i, k) = zfice(i, k)*pqlwp(i, k)
267
      ENDDO
268
    ENDDO
269
  ENDIF
270
271
72
  IF (ok_cdnc) THEN
272
273
    ! --we compute cloud properties as a function of the aerosol load
274
275
    DO k = 1, klev
276
      DO i = 1, klon
277
        ! Formula "D" of Boucher and Lohmann, Tellus, 1995
278
        ! Cloud droplet number concentration (CDNC) is restricted
279
        ! to be within [20, 1000 cm^3]
280
281
        ! --pre-industrial case
282
        cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
283
          1.E-4))/log(10.))*1.E6 !-m-3
284
        cdnc_pi(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc_pi(i,k)))
285
286
      ENDDO
287
    ENDDO
288
289
    !--flag_aerosol=7 => MACv2SP climatology
290
    !--in this case there is an enhancement factor
291
    IF (flag_aerosol .EQ. 7) THEN
292
293
      !--present-day
294
      DO k = 1, klev
295
        DO i = 1, klon
296
          cdnc(i, k) = cdnc_pi(i,k)*dNovrN(i)
297
        ENDDO
298
      ENDDO
299
300
    !--standard case
301
    ELSE
302
303
      DO k = 1, klev
304
        DO i = 1, klon
305
306
          ! Formula "D" of Boucher and Lohmann, Tellus, 1995
307
          ! Cloud droplet number concentration (CDNC) is restricted
308
          ! to be within [20, 1000 cm^3]
309
310
          ! --present-day case
311
          cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
312
            1.E-4))/log(10.))*1.E6 !-m-3
313
          cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k)))
314
315
        ENDDO
316
      ENDDO
317
318
    ENDIF !--flag_aerosol
319
320
    !--computing cloud droplet size
321
    DO k = 1, klev
322
      DO i = 1, klon
323
324
        ! --present-day case
325
        rad_chaud(i, k) = 1.1*((pqlwp(i,k)*pplay(i, &
326
          k)/(rd*t(i,k)))/(4./3*rpi*1000.*cdnc(i,k)))**(1./3.)
327
        rad_chaud(i, k) = max(rad_chaud(i,k)*1.E6, 5.)
328
329
        ! --pre-industrial case
330
        rad_chaud_pi(i, k) = 1.1*((pqlwp(i,k)*pplay(i, &
331
          k)/(rd*t(i,k)))/(4./3.*rpi*1000.*cdnc_pi(i,k)))**(1./3.)
332
        rad_chaud_pi(i, k) = max(rad_chaud_pi(i,k)*1.E6, 5.)
333
334
        ! --pre-industrial case
335
        ! --liquid/ice cloud water paths:
336
        IF (pclc(i,k)<=seuil_neb) THEN
337
338
          pcldtaupi(i, k) = 0.0
339
340
        ELSE
341
342
          zflwp_var = 1000.*(1.-zfice(i,k))*pqlwp(i, k)/pclc(i, k)* &
343
            rhodz(i, k)
344
          zfiwp_var = 1000.*zfice(i, k)*pqlwp(i, k)/pclc(i, k)*rhodz(i, k)
345
          ! Calculation of ice cloud effective radius in micron
346
          IF (iflag_rei .EQ. 1) THEN
347
            ! when we account for precipitation in the radiation scheme,
348
            ! It is recommended to use the rei formula from Sun and Rikkus 1999 with a revision
349
            ! from Sun 2001 (as in the IFS model)
350
            iwc=zfice(i, k)*pqlwp(i, k)/pclc(i,k)*zrho(i,k)*1000. !in cloud ice water content in g/m3
351
            dei=(1.2351+0.0105*(t(i,k)-273.15))*(45.8966*(iwc**0.2214) + &
352
               & 0.7957*(iwc**0.2535)*(t(i,k)-83.15))
353
            !deimax=155.0
354
            !deimin=20.+40*cos(abs(latitude_deg(i))/180.*RPI)
355
            !Etienne: deimax and deimin controled by rei_max and rei_min in physiq.def
356
            deimax=rei_max*2.0
357
            deimin=2.0*rei_min+40*cos(abs(latitude_deg(i))/180.*RPI)
358
            dei=min(dei,deimax)
359
            dei=max(dei,deimin)
360
            rei=3.*sqrt(3.)/8.*dei
361
           ELSE
362
            ! Default
363
            ! for ice clouds: as a function of the ambiant temperature
364
            ! [formula used by Iacobellis and Somerville (2000), with an
365
            ! asymptotical value of 3.5 microns at T<-81.4 C added to be
366
            ! consistent with observations of Heymsfield et al. 1986]:
367
            ! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
368
            ! rei_max=61.29
369
            tc = t(i, k) - 273.15
370
            rei = d_rei_dt*tc + rei_max
371
            IF (tc<=-81.4) rei = rei_min
372
           ENDIF
373
374
          ! -- cloud optical thickness :
375
          ! [for liquid clouds, traditional formula,
376
          ! for ice clouds, Ebert & Curry (1992)]
377
378
          IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
379
          pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
380
            zfiwp_var*(3.448E-03+2.431/rei)
381
382
        ENDIF
383
384
      ENDDO
385
    ENDDO
386
387
  ELSE !--not ok_cdnc
388
389
    ! -prescribed cloud droplet radius
390
391
288
    DO k = 1, min(3, klev)
392
214992
      DO i = 1, klon
393
214704
        rad_chaud(i, k) = rad_chau2
394
214920
        rad_chaud_pi(i, k) = rad_chau2
395
      ENDDO
396
    ENDDO
397
2664
    DO k = min(3, klev) + 1, klev
398
2579112
      DO i = 1, klon
399
2576448
        rad_chaud(i, k) = rad_chau1
400
2579040
        rad_chaud_pi(i, k) = rad_chau1
401
      ENDDO
402
    ENDDO
403
404
  ENDIF !--ok_cdnc
405
406
  ! --computation of cloud optical depth and emissivity
407
  ! --in the general case
408
409
2880
  DO k = 1, klev
410
2794032
    DO i = 1, klon
411
412
2791152
      IF (pclc(i,k)<=seuil_neb) THEN
413
414
        ! effective cloud droplet radius (microns) for liquid water clouds:
415
        ! For output diagnostics cloud droplet effective radius [um]
416
        ! we multiply here with f * xl (fraction of liquid water
417
        ! clouds in the grid cell) to avoid problems in the averaging of the
418
        ! output.
419
        ! In the output of IOIPSL, derive the REAL cloud droplet
420
        ! effective radius as re/fl
421
422
2170851
        fl(i, k) = seuil_neb*(1.-zfice(i,k))
423
2170851
        re(i, k) = rad_chaud(i, k)*fl(i, k)
424
        rel = 0.
425
        rei = 0.
426
2170851
        pclc(i, k) = 0.0
427
2170851
        pcltau(i, k) = 0.0
428
2170851
        pclemi(i, k) = 0.0
429
430
      ELSE
431
432
        ! -- liquid/ice cloud water paths:
433
434
620301
        zflwp_var = 1000.*(1.-zfice(i,k))*pqlwp(i, k)/pclc(i, k)*rhodz(i, k)
435
620301
        zfiwp_var = 1000.*zfice(i, k)*pqlwp(i, k)/pclc(i, k)*rhodz(i, k)
436
437
        ! effective cloud droplet radius (microns) for liquid water clouds:
438
        ! For output diagnostics cloud droplet effective radius [um]
439
        ! we multiply here with f * xl (fraction of liquid water
440
        ! clouds in the grid cell) to avoid problems in the averaging of the
441
        ! output.
442
        ! In the output of IOIPSL, derive the REAL cloud droplet
443
        ! effective radius as re/fl
444
445
620301
        fl(i, k) = pclc(i, k)*(1.-zfice(i,k))
446
620301
        re(i, k) = rad_chaud(i, k)*fl(i, k)
447
448
        rel = rad_chaud(i, k)
449
450
        ! Calculation of ice cloud effective radius in micron
451
452
453
620301
        IF (iflag_rei .GT. 0) THEN
454
455
            ! when we account for precipitation in the radiation scheme,
456
            ! we use the rei formula from Sun and Rikkus 1999 with a revision
457
            ! from Sun 2001 (as in the IFS model)
458
            iwc=zfice(i, k)*pqlwp(i, k)/pclc(i,k)*zrho(i,k)*1000. !in cloud ice water content in g/m3
459
            dei=(1.2351+0.0105*(t(i,k)-273.15))*(45.8966*(iwc**0.2214) + &
460
               &0.7957*(iwc**0.2535)*(t(i,k)-83.15))
461
            !deimax=155.0
462
            !deimin=20.+40*cos(abs(latitude_deg(i))/180.*RPI)
463
            !Etienne: deimax and deimin controled by rei_max and rei_min in physiq.def
464
            deimax=rei_max*2.0
465
            deimin=2.0*rei_min+40*cos(abs(latitude_deg(i))/180.*RPI)
466
            dei=min(dei,deimax)
467
            dei=max(dei,deimin)
468
            rei=3.*sqrt(3.)/8.*dei
469
470
        ELSE
471
            ! Default
472
            ! for ice clouds: as a function of the ambiant temperature
473
            ! [formula used by Iacobellis and Somerville (2000), with an
474
            ! asymptotical value of 3.5 microns at T<-81.4 C added to be
475
            ! consistent with observations of Heymsfield et al. 1986]:
476
            ! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
477
            ! rei_max=61.29
478
620301
            tc = t(i, k) - 273.15
479
620301
            rei = d_rei_dt*tc + rei_max
480
620301
            IF (tc<=-81.4) rei = rei_min
481
        ENDIF
482
        ! -- cloud optical thickness :
483
        ! [for liquid clouds, traditional formula,
484
        ! for ice clouds, Ebert & Curry (1992)]
485
486
620301
        IF (zflwp_var==0.) rel = 1.
487

620301
        IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
488
        pcltau(i, k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/ &
489
620301
          rei)
490
491
        ! -- cloud infrared emissivity:
492
        ! [the broadband infrared absorption coefficient is PARAMETERized
493
        ! as a function of the effective cld droplet radius]
494
        ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
495
496
620301
        k_ice = k_ice0 + 1.0/rei
497
498
620301
        pclemi(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
499
500
      ENDIF
501
502
2791152
      reice(i, k) = rei
503
504
2791152
      xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
505
2793960
      xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
506
507
    ENDDO
508
  ENDDO
509
510
  ! --if cloud droplet radius is fixed, then pcldtaupi=pcltau
511
512
72
  IF (.NOT. ok_cdnc) THEN
513
2880
    DO k = 1, klev
514
2794032
      DO i = 1, klon
515
2791152
        pcldtaupi(i, k) = pcltau(i, k)
516
2793960
        reice_pi(i, k) = reice(i, k)
517
      ENDDO
518
    ENDDO
519
  ENDIF
520
521
2880
  DO k = 1, klev
522
2794032
    DO i = 1, klon
523
2791152
      reliq(i, k) = rad_chaud(i, k)
524
2791152
      reliq_pi(i, k) = rad_chaud_pi(i, k)
525
2793960
      reice_pi(i, k) = reice(i, k)
526
    ENDDO
527
  ENDDO
528
529
  ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
530
  ! IM cf. CR:test: calcul prenant ou non en compte le recouvrement
531
  ! initialisations
532
533
71640
  DO i = 1, klon
534
71568
    zclear(i) = 1.
535
71568
    zcloud(i) = 0.
536
71568
    zcloudh(i) = 0.
537
71568
    zcloudm(i) = 0.
538
71568
    zcloudl(i) = 0.
539
71568
    pch(i) = 1.0
540
71568
    pcm(i) = 1.0
541
71568
    pcl(i) = 1.0
542
71640
    pctlwp(i) = 0.0
543
  ENDDO
544
545
  ! --calculation of liquid water path
546
547
2880
  DO k = klev, 1, -1
548
2794032
    DO i = 1, klon
549
2793960
      pctlwp(i) = pctlwp(i) + pqlwp(i, k)*rhodz(i, k)
550
    ENDDO
551
  ENDDO
552
553
  ! --calculation of cloud properties with cloud overlap
554
555
  IF (novlp==1) THEN
556
2880
    DO k = klev, 1, -1
557
2794032
      DO i = 1, klon
558
        zclear(i) = zclear(i)*(1.-max(pclc(i,k),zcloud(i)))/(1.-min(real( &
559
2791152
          zcloud(i),kind=8),1.-zepsec))
560
2791152
        pct(i) = 1. - zclear(i)
561
2791152
        IF (paprs(i,k)<prmhc) THEN
562
          pch(i) = pch(i)*(1.-max(pclc(i,k),zcloudh(i)))/(1.-min(real(zcloudh &
563
1869880
            (i),kind=8),1.-zepsec))
564
1869880
          zcloudh(i) = pclc(i, k)
565

921272
        ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
566
          pcm(i) = pcm(i)*(1.-max(pclc(i,k),zcloudm(i)))/(1.-min(real(zcloudm &
567
247328
            (i),kind=8),1.-zepsec))
568
247328
          zcloudm(i) = pclc(i, k)
569
673944
        ELSE IF (paprs(i,k)>=prlmc) THEN
570
          pcl(i) = pcl(i)*(1.-max(pclc(i,k),zcloudl(i)))/(1.-min(real(zcloudl &
571
673944
            (i),kind=8),1.-zepsec))
572
673944
          zcloudl(i) = pclc(i, k)
573
        ENDIF
574
2793960
        zcloud(i) = pclc(i, k)
575
      ENDDO
576
    ENDDO
577
  ELSE IF (novlp==2) THEN
578
    DO k = klev, 1, -1
579
      DO i = 1, klon
580
        zcloud(i) = max(pclc(i,k), zcloud(i))
581
        pct(i) = zcloud(i)
582
        IF (paprs(i,k)<prmhc) THEN
583
          pch(i) = min(pclc(i,k), pch(i))
584
        ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
585
          pcm(i) = min(pclc(i,k), pcm(i))
586
        ELSE IF (paprs(i,k)>=prlmc) THEN
587
          pcl(i) = min(pclc(i,k), pcl(i))
588
        ENDIF
589
      ENDDO
590
    ENDDO
591
  ELSE IF (novlp==3) THEN
592
    DO k = klev, 1, -1
593
      DO i = 1, klon
594
        zclear(i) = zclear(i)*(1.-pclc(i,k))
595
        pct(i) = 1 - zclear(i)
596
        IF (paprs(i,k)<prmhc) THEN
597
          pch(i) = pch(i)*(1.0-pclc(i,k))
598
        ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
599
          pcm(i) = pcm(i)*(1.0-pclc(i,k))
600
        ELSE IF (paprs(i,k)>=prlmc) THEN
601
          pcl(i) = pcl(i)*(1.0-pclc(i,k))
602
        ENDIF
603
      ENDDO
604
    ENDDO
605
  ENDIF
606
607
71640
  DO i = 1, klon
608
71568
    pch(i) = 1. - pch(i)
609
71568
    pcm(i) = 1. - pcm(i)
610
71640
    pcl(i) = 1. - pcl(i)
611
  ENDDO
612
613
  ! ========================================================
614
  ! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
615
  ! ========================================================
616
  ! change by Nicolas Yan (LSCE)
617
  ! Cloud Droplet Number Concentration (CDNC) : 3D variable
618
  ! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
619
  ! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
620
  ! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
621
  ! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
622
623
72
  IF (ok_cdnc) THEN
624
625
    DO k = 1, klev
626
      DO i = 1, klon
627
        phase3d(i, k) = 1 - zfice(i, k)
628
        IF (pclc(i,k)<=seuil_neb) THEN
629
          lcc3d(i, k) = seuil_neb*phase3d(i, k)
630
        ELSE
631
          lcc3d(i, k) = pclc(i, k)*phase3d(i, k)
632
        ENDIF
633
        scdnc(i, k) = lcc3d(i, k)*cdnc(i, k) ! m-3
634
      ENDDO
635
    ENDDO
636
637
    DO i = 1, klon
638
      lcc(i) = 0.
639
      reffclwtop(i) = 0.
640
      cldncl(i) = 0.
641
      IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1.
642
      IF (novlp.EQ.2) tcc(i) = 0.
643
    ENDDO
644
645
    DO i = 1, klon
646
      DO k = klev - 1, 1, -1 !From TOA down
647
648
          ! Test, if the cloud optical depth exceeds the necessary
649
          ! threshold:
650
651
        IF (pcltau(i,k)>thres_tau .AND. pclc(i,k)>thres_neb) THEN
652
653
          IF (novlp.EQ.2) THEN
654
            IF (first) THEN
655
              WRITE (*, *) 'Hypothese de recouvrement: MAXIMUM'
656
              first = .FALSE.
657
            ENDIF
658
            flag_max = -1.
659
            ftmp(i) = max(tcc(i), pclc(i,k))
660
          ENDIF
661
662
          IF (novlp.EQ.3) THEN
663
            IF (first) THEN
664
              WRITE (*, *) 'Hypothese de recouvrement: RANDOM'
665
              first = .FALSE.
666
            ENDIF
667
            flag_max = 1.
668
            ftmp(i) = tcc(i)*(1-pclc(i,k))
669
          ENDIF
670
671
          IF (novlp.EQ.1) THEN
672
            IF (first) THEN
673
              WRITE (*, *) 'Hypothese de recouvrement: MAXIMUM_ &
674
                &                                             &
675
                &                                          RANDOM'
676
              first = .FALSE.
677
            ENDIF
678
            flag_max = 1.
679
            ftmp(i) = tcc(i)*(1.-max(pclc(i,k),pclc(i,k+1)))/(1.-min(pclc(i, &
680
              k+1),1.-thres_neb))
681
          ENDIF
682
          ! Effective radius of cloud droplet at top of cloud (m)
683
          reffclwtop(i) = reffclwtop(i) + rad_chaud(i, k)*1.0E-06*phase3d(i, &
684
            k)*(tcc(i)-ftmp(i))*flag_max
685
          ! CDNC at top of cloud (m-3)
686
          cldncl(i) = cldncl(i) + cdnc(i, k)*phase3d(i, k)*(tcc(i)-ftmp(i))* &
687
            flag_max
688
          ! Liquid Cloud Content at top of cloud
689
          lcc(i) = lcc(i) + phase3d(i, k)*(tcc(i)-ftmp(i))*flag_max
690
          ! Total Cloud Content at top of cloud
691
          tcc(i) = ftmp(i)
692
693
        ENDIF ! is there a visible, not-too-small cloud?
694
      ENDDO ! loop over k
695
696
      IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1. - tcc(i)
697
698
    ENDDO ! loop over i
699
700
    ! ! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC
701
    ! REFFCLWS)
702
    DO i = 1, klon
703
      DO k = 1, klev
704
        ! Weight to be used for outputs: eau_liquide*couverture nuageuse
705
        lcc3dcon(i, k) = rnebcon(i, k)*phase3d(i, k)*clwcon(i, k) ! eau liquide convective
706
        lcc3dstra(i, k) = pclc(i, k)*pqlwp(i, k)*phase3d(i, k)
707
        lcc3dstra(i, k) = lcc3dstra(i, k) - lcc3dcon(i, k) ! eau liquide stratiforme
708
        lcc3dstra(i, k) = max(lcc3dstra(i,k), 0.0)
709
        !FC pour la glace (CAUSES)
710
        icc3dcon(i, k) = rnebcon(i, k)*(1-phase3d(i, k))*clwcon(i, k) !  glace convective
711
        icc3dstra(i, k)= pclc(i, k)*pqlwp(i, k)*(1-phase3d(i, k))
712
        icc3dstra(i, k) = icc3dstra(i, k) - icc3dcon(i, k) ! glace stratiforme
713
        icc3dstra(i, k) = max( icc3dstra(i, k), 0.0)
714
        !FC (CAUSES)
715
716
        ! Compute cloud droplet radius as above in meter
717
        radius = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3*rpi*1000.* &
718
          cdnc(i,k)))**(1./3.)
719
        radius = max(radius, 5.E-6)
720
        ! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
721
        reffclwc(i, k) = radius
722
        reffclwc(i, k) = reffclwc(i, k)*lcc3dcon(i, k)
723
        ! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
724
        reffclws(i, k) = radius
725
        reffclws(i, k) = reffclws(i, k)*lcc3dstra(i, k)
726
      ENDDO !klev
727
    ENDDO !klon
728
729
    ! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
730
731
    DO i = 1, klon
732
      cldnvi(i) = 0.
733
      lcc_integrat(i) = 0.
734
      height(i) = 0.
735
      DO k = 1, klev
736
        cldnvi(i) = cldnvi(i) + cdnc(i, k)*lcc3d(i, k)*dh(i, k)
737
        lcc_integrat(i) = lcc_integrat(i) + lcc3d(i, k)*dh(i, k)
738
        height(i) = height(i) + dh(i, k)
739
      ENDDO ! klev
740
      lcc_integrat(i) = lcc_integrat(i)/height(i)
741
      IF (lcc_integrat(i)<=1.0E-03) THEN
742
        cldnvi(i) = cldnvi(i)*lcc(i)/seuil_neb
743
      ELSE
744
        cldnvi(i) = cldnvi(i)*lcc(i)/lcc_integrat(i)
745
      ENDIF
746
    ENDDO ! klon
747
748
    DO i = 1, klon
749
      DO k = 1, klev
750
        IF (scdnc(i,k)<=0.0) scdnc(i, k) = 0.0
751
        IF (reffclws(i,k)<=0.0) reffclws(i, k) = 0.0
752
        IF (reffclwc(i,k)<=0.0) reffclwc(i, k) = 0.0
753
        IF (lcc3d(i,k)<=0.0) lcc3d(i, k) = 0.0
754
        IF (lcc3dcon(i,k)<=0.0) lcc3dcon(i, k) = 0.0
755
        IF (lcc3dstra(i,k)<=0.0) lcc3dstra(i, k) = 0.0
756
!FC (CAUSES)
757
        IF (icc3dcon(i,k)<=0.0) icc3dcon(i, k) = 0.0
758
        IF (icc3dstra(i,k)<=0.0) icc3dstra(i, k) = 0.0
759
!FC (CAUSES)
760
      ENDDO
761
      IF (reffclwtop(i)<=0.0) reffclwtop(i) = 0.0
762
      IF (cldncl(i)<=0.0) cldncl(i) = 0.0
763
      IF (cldnvi(i)<=0.0) cldnvi(i) = 0.0
764
      IF (lcc(i)<=0.0) lcc(i) = 0.0
765
    ENDDO
766
767
  ENDIF !ok_cdnc
768
769
72
  first=.false. !to be sure
770
771
72
  RETURN
772
773
END SUBROUTINE newmicro