source: LMDZ6/branches/Amaury_dev/libf/phylmd/nuage.F90 @ 5151

Last change on this file since 5151 was 5144, checked in by abarral, 7 weeks ago

Put YOMCST.h into modules

  • 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: 14.0 KB
Line 
1! $Id: nuage.F90 5144 2024-07-29 21:01:04Z abarral $
2
3SUBROUTINE nuage(paprs, pplay, t, pqlwp, picefra, pclc, pcltau, pclemi, pch, pcl, pcm, &
4        pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, distcltop, &
5        temp_cltop, cldtaupi, re, fl)
6  USE dimphy
7  USE lmdz_lscp_tools, ONLY: icefrac_lscp
8  USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
9  USE lmdz_lscp_ini, ONLY: iflag_t_glace
10  USE phys_local_var_mod, ONLY: ptconv
11  USE lmdz_clesphys
12  USE lmdz_nuage_params ! JBM 3/14
13  USE lmdz_yomcst
14
15  IMPLICIT NONE
16  ! ======================================================================
17  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
18  ! Objet: Calculer epaisseur optique et emmissivite des nuages
19  ! ======================================================================
20  ! Arguments:
21  ! t-------input-R-temperature
22  ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
23  ! picefra--inout-R-fraction de glace dans les nuages (-)
24  ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
25  ! ok_aie--input-L-apply aerosol indirect effect or not
26  ! mass_solu_aero-----input-R-total mass concentration for all soluble
27  ! aerosols[ug/m^3]
28  ! mass_solu_aero_pi--input-R-dito, pre-industrial value
29  ! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
30  ! bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
31
32  ! cldtaupi-output-R-pre-industrial value of cloud optical thickness,
33  ! needed for the diagnostics of the aerosol indirect
34  ! radiative forcing (see radlwsw)
35  ! re------output-R-Cloud droplet effective radius multiplied by fl [um]
36  ! fl------output-R-Denominator to re, introduced to avoid problems in
37  ! the averaging of the output. fl is the fraction of liquid
38  ! water clouds within a grid cell
39
40  ! pcltau--output-R-epaisseur optique des nuages
41  ! pclemi--output-R-emissivite des nuages (0 a 1)
42  ! ======================================================================
43
44  REAL paprs(klon, klev + 1), pplay(klon, klev)
45  REAL t(klon, klev)
46
47  REAL pclc(klon, klev)
48  REAL pqlwp(klon, klev), picefra(klon, klev)
49  REAL pcltau(klon, klev), pclemi(klon, klev)
50
51  REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
52  REAL distcltop(klon, klev)
53  REAL temp_cltop(klon, klev)
54  LOGICAL lo
55
56  REAL cetahb, cetamb
57  PARAMETER (cetahb = 0.45, cetamb = 0.80)
58
59  INTEGER i, k
60  REAL zflwp, zradef, zfice(klon), zmsac
61
62  REAL radius, rad_chaud
63  ! JBM (3/14) parameters already defined in nuage.h:
64  ! REAL rad_froid, rad_chau1, rad_chau2
65  ! PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
66  ! cc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
67  ! sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
68  REAL coef, coef_froi, coef_chau
69  PARAMETER (coef_chau = 0.13, coef_froi = 0.09)
70  REAL seuil_neb
71  PARAMETER (seuil_neb = 0.001)
72  ! JBM (3/14) nexpo is replaced by exposant_glace
73  ! REAL nexpo ! exponentiel pour glace/eau
74  ! PARAMETER (nexpo=6.)
75  REAL, PARAMETER :: t_glace_min_old = 258.
76  INTEGER, PARAMETER :: exposant_glace_old = 6
77
78
79  ! jq for the aerosol indirect effect
80  ! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
81  ! jq
82  LOGICAL ok_aie ! Apply AIE or not?
83
84  REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols[ug m-3]
85  REAL mass_solu_aero_pi(klon, klev) ! - " - pre-industrial value
86  REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
87  REAL re(klon, klev) ! cloud droplet effective radius [um]
88  REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
89  REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
90
91  REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds
92  ! within the grid cell)
93
94  REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
95
96  REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
97  REAl dzfice(klon)
98  ! jq-end
99
100  ! cc      PARAMETER (nexpo=1)
101
102  ! Calculer l'epaisseur optique et l'emmissivite des nuages
103
104  DO k = 1, klev
105    IF (iflag_t_glace==0) THEN
106      DO i = 1, klon
107        zfice(i) = 1.0 - (t(i, k) - t_glace_min_old) / (273.13 - t_glace_min_old)
108        zfice(i) = min(max(zfice(i), 0.0), 1.0)
109        zfice(i) = zfice(i)**exposant_glace_old
110      ENDDO
111    ELSE ! of IF (iflag_t_glace.EQ.0)
112      ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
113      !       zfice(i) = icefrac_lsc(t(i,k), t_glace_min, &
114      !                           t_glace_max, exposant_glace)
115      IF (ok_new_lscp) THEN
116        CALL icefrac_lscp(klon, t(:, k), iflag_ice_thermo, distcltop(:, k), temp_cltop(:, k), zfice(:), dzfice(:))
117      ELSE
118        CALL icefrac_lsc(klon, t(:, k), pplay(:, k) / paprs(:, 1), zfice(:))
119
120      ENDIF
121
122      IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
123        ! EV: take the ice fraction directly from the lscp code
124        ! consistent only for non convective grid points
125        ! critical for mixed phase clouds
126        DO i = 1, klon
127          IF (.NOT. ptconv(i, k)) THEN
128            zfice(i) = picefra(i, k)
129          ENDIF
130        ENDDO
131      ENDIF
132
133    ENDIF
134
135    DO i = 1, klon
136      rad_chaud = rad_chau1
137      IF (k<=3) rad_chaud = rad_chau2
138
139      pclc(i, k) = max(pclc(i, k), seuil_neb)
140      zflwp = 1000. * pqlwp(i, k) / rg / pclc(i, k) * (paprs(i, k) - paprs(i, k + 1))
141
142      IF (ok_aie) THEN
143        ! Formula "D" of Boucher and Lohmann, Tellus, 1995
144
145        cdnc(i, k) = 10.**(bl95_b0 + bl95_b1 * log(max(mass_solu_aero(i, k), &
146                1.E-4)) / log(10.)) * 1.E6 !-m-3
147        ! Cloud droplet number concentration (CDNC) is restricted
148        ! to be within [20, 1000 cm^3]
149
150        cdnc(i, k) = min(1000.E6, max(20.E6, cdnc(i, k)))
151        cdnc_pi(i, k) = 10.**(bl95_b0 + bl95_b1 * log(max(mass_solu_aero_pi(i, k), &
152                1.E-4)) / log(10.)) * 1.E6 !-m-3
153        cdnc_pi(i, k) = min(1000.E6, max(20.E6, cdnc_pi(i, k)))
154
155
156        ! air density: pplay(i,k) / (RD * zT(i,k))
157        ! factor 1.1: derive effective radius from volume-mean radius
158        ! factor 1000 is the water density
159        ! _chaud means that this is the CDR for liquid water clouds
160
161        rad_chaud = 1.1 * ((pqlwp(i, k) * pplay(i, k) / (rd * t(i, k))) / (4. / 3. * rpi * 1000. &
162                * cdnc(i, k)))**(1. / 3.)
163
164        ! Convert to um. CDR shall be at least 3 um.
165
166        rad_chaud = max(rad_chaud * 1.E6, 3.)
167
168        ! For output diagnostics
169
170        ! Cloud droplet effective radius [um]
171
172        ! we multiply here with f * xl (fraction of liquid water
173        ! clouds in the grid cell) to avoid problems in the
174        ! averaging of the output.
175        ! In the output of IOIPSL, derive the real cloud droplet
176        ! effective radius as re/fl
177
178        fl(i, k) = pclc(i, k) * (1. - zfice(i))
179        re(i, k) = rad_chaud * fl(i, k)
180
181        ! Pre-industrial cloud opt thickness
182
183        ! "radius" is calculated as rad_chaud above (plus the
184        ! ice cloud contribution) but using cdnc_pi instead of
185        ! cdnc.
186        radius = max(1.1E6 * ((pqlwp(i, k) * pplay(i, k) / (rd * t(i, k))) / (4. / 3. * rpi * &
187                1000. * cdnc_pi(i, k)))**(1. / 3.), 3.) * (1. - zfice(i)) + rad_froid * zfice(i)
188        cldtaupi(i, k) = 3.0 / 2.0 * zflwp / radius
189      END IF ! ok_aie
190
191      radius = rad_chaud * (1. - zfice(i)) + rad_froid * zfice(i)
192      coef = coef_chau * (1. - zfice(i)) + coef_froi * zfice(i)
193      pcltau(i, k) = 3.0 / 2.0 * zflwp / radius
194      pclemi(i, k) = 1.0 - exp(-coef * zflwp)
195      lo = (pclc(i, k)<=seuil_neb)
196      IF (lo) pclc(i, k) = 0.0
197      IF (lo) pcltau(i, k) = 0.0
198      IF (lo) pclemi(i, k) = 0.0
199
200      IF (.NOT. ok_aie) cldtaupi(i, k) = pcltau(i, k)
201    END DO
202  END DO
203  ! cc      DO k = 1, klev
204  ! cc      DO i = 1, klon
205  ! cc         t(i,k) = t(i,k)
206  ! cc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
207  ! cc         lo = pclc(i,k) .GT. (2.*1.e-5)
208  ! cc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
209  ! cc     .          /(rg*pclc(i,k))
210  ! cc         zradef = 10.0 + (1.-sigs(k))*45.0
211  ! cc         pcltau(i,k) = 1.5 * zflwp / zradef
212  ! cc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
213  ! cc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
214  ! cc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
215  ! cc         if (.NOT.lo) pclc(i,k) = 0.0
216  ! cc         if (.NOT.lo) pcltau(i,k) = 0.0
217  ! cc         if (.NOT.lo) pclemi(i,k) = 0.0
218  ! cc      ENDDO
219  ! cc      ENDDO
220  ! ccccc      PRINT*, 'pas de nuage dans le rayonnement'
221  ! ccccc      DO k = 1, klev
222  ! ccccc      DO i = 1, klon
223  ! ccccc         pclc(i,k) = 0.0
224  ! ccccc         pcltau(i,k) = 0.0
225  ! ccccc         pclemi(i,k) = 0.0
226  ! ccccc      ENDDO
227  ! ccccc      ENDDO
228
229  ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
230
231  DO i = 1, klon
232    pct(i) = 1.0
233    pch(i) = 1.0
234    pcm(i) = 1.0
235    pcl(i) = 1.0
236    pctlwp(i) = 0.0
237  END DO
238
239  DO k = klev, 1, -1
240    DO i = 1, klon
241      pctlwp(i) = pctlwp(i) + pqlwp(i, k) * (paprs(i, k) - paprs(i, k + 1)) / rg
242      pct(i) = pct(i) * (1.0 - pclc(i, k))
243      IF (pplay(i, k)<=cetahb * paprs(i, 1)) pch(i) = pch(i) * (1.0 - pclc(i, k))
244      IF (pplay(i, k)>cetahb * paprs(i, 1) .AND. pplay(i, k)<=cetamb * paprs(i, 1)) &
245              pcm(i) = pcm(i) * (1.0 - pclc(i, k))
246      IF (pplay(i, k)>cetamb * paprs(i, 1)) pcl(i) = pcl(i) * (1.0 - pclc(i, k))
247    END DO
248  END DO
249
250  DO i = 1, klon
251    pct(i) = 1. - pct(i)
252    pch(i) = 1. - pch(i)
253    pcm(i) = 1. - pcm(i)
254    pcl(i) = 1. - pcl(i)
255  END DO
256
257END SUBROUTINE nuage
258SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq)
259  USE dimphy
260  USE lmdz_yomcst
261
262  IMPLICIT NONE
263
264  ! Laurent Li (LMD/CNRS), le 12 octobre 1998
265  ! (adaptation du code ECMWF)
266
267  ! Dans certains cas, le schema pronostique des nuages n'est
268  ! pas suffisament performant. On a donc besoin de diagnostiquer
269  ! ces nuages. Je dois avouer que c'est une frustration.
270
271  ! Arguments d'entree:
272  REAL paprs(klon, klev + 1) ! pression (Pa) a inter-couche
273  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
274  REAL t(klon, klev) ! temperature (K)
275  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
276  REAL rain(klon) ! pluie convective (kg/m2/s)
277  REAL snow(klon) ! neige convective (kg/m2/s)
278  INTEGER ktop(klon) ! sommet de la convection
279  INTEGER kbot(klon) ! bas de la convection
280
281  ! Arguments de sortie:
282  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
283  REAL dialiq(klon, klev) ! eau liquide nuageuse
284
285  ! Constantes ajustables:
286  REAL canva, canvb, canvh
287  PARAMETER (canva = 2.0, canvb = 0.3, canvh = 0.4)
288  REAL cca, ccb, ccc
289  PARAMETER (cca = 0.125, ccb = 1.5, ccc = 0.8)
290  REAL ccfct, ccscal
291  PARAMETER (ccfct = 0.400)
292  PARAMETER (ccscal = 1.0E+11)
293  REAL cetahb, cetamb
294  PARAMETER (cetahb = 0.45, cetamb = 0.80)
295  REAL cclwmr
296  PARAMETER (cclwmr = 1.E-04)
297  REAL zepscr
298  PARAMETER (zepscr = 1.0E-10)
299
300  ! Variables locales:
301  INTEGER i, k
302  REAL zcc(klon)
303
304  ! Initialisation:
305
306  DO k = 1, klev
307    DO i = 1, klon
308      diafra(i, k) = 0.0
309      dialiq(i, k) = 0.0
310    END DO
311  END DO
312
313  DO i = 1, klon ! Calculer la fraction nuageuse
314    zcc(i) = 0.0
315    IF ((rain(i) + snow(i))>0.) THEN
316      zcc(i) = cca * log(max(zepscr, (rain(i) + snow(i)) * ccscal)) - ccb
317      zcc(i) = min(ccc, max(0.0, zcc(i)))
318    END IF
319  END DO
320
321  DO i = 1, klon ! pour traiter les enclumes
322    diafra(i, ktop(i)) = max(diafra(i, ktop(i)), zcc(i) * ccfct)
323    IF ((zcc(i)>=canvh) .AND. (pplay(i, ktop(i))<=cetahb * paprs(i, &
324            1))) diafra(i, ktop(i)) = max(diafra(i, ktop(i)), max(zcc(&
325            i) * ccfct, canva * (zcc(i) - canvb)))
326    dialiq(i, ktop(i)) = cclwmr * diafra(i, ktop(i))
327  END DO
328
329  DO k = 1, klev ! nuages convectifs (sauf enclumes)
330    DO i = 1, klon
331      IF (k<ktop(i) .AND. k>=kbot(i)) THEN
332        diafra(i, k) = max(diafra(i, k), zcc(i) * ccfct)
333        dialiq(i, k) = cclwmr * diafra(i, k)
334      END IF
335    END DO
336  END DO
337
338END SUBROUTINE diagcld1
339SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
340  USE dimphy
341  USE lmdz_yoethf
342  USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
343  USE lmdz_yomcst
344
345  IMPLICIT NONE
346
347  ! Arguments d'entree:
348  REAL paprs(klon, klev + 1) ! pression (Pa) a inter-couche
349  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
350  REAL t(klon, klev) ! temperature (K)
351  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
352
353  ! Arguments de sortie:
354  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
355  REAL dialiq(klon, klev) ! eau liquide nuageuse
356
357  REAL cetamb
358  PARAMETER (cetamb = 0.80)
359  REAL cloia, cloib, cloic, cloid
360  PARAMETER (cloia = 1.0E+02, cloib = -10.00, cloic = -0.6, cloid = 5.0)
361  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
362  REAL rgammas
363  PARAMETER (rgammas = 0.05)
364  REAL crhl
365  PARAMETER (crhl = 0.15)
366  ! cc      PARAMETER (CRHL=0.70)
367  REAL t_coup
368  PARAMETER (t_coup = 234.0)
369
370  ! Variables locales:
371  INTEGER i, k, kb, invb(klon)
372  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
373  REAL zdelta, zcor
374
375  ! Initialisation:
376
377  DO k = 1, klev
378    DO i = 1, klon
379      diafra(i, k) = 0.0
380      dialiq(i, k) = 0.0
381    END DO
382  END DO
383
384  DO i = 1, klon
385    invb(i) = klev
386    zdthmin(i) = 0.0
387  END DO
388
389  DO k = 2, klev / 2 - 1
390    DO i = 1, klon
391      zdthdp = (t(i, k) - t(i, k + 1)) / (pplay(i, k) - pplay(i, k + 1)) - &
392              rd * 0.5 * (t(i, k) + t(i, k + 1)) / rcpd / paprs(i, k + 1)
393      zdthdp = zdthdp * cloia
394      IF (pplay(i, k)>cetamb * paprs(i, 1) .AND. zdthdp<zdthmin(i)) THEN
395        zdthmin(i) = zdthdp
396        invb(i) = k
397      END IF
398    END DO
399  END DO
400
401  DO i = 1, klon
402    kb = invb(i)
403    IF (thermcep) THEN
404      zdelta = max(0., sign(1., rtt - t(i, kb)))
405      zqs = r2es * foeew(t(i, kb), zdelta) / pplay(i, kb)
406      zqs = min(0.5, zqs)
407      zcor = 1. / (1. - retv * zqs)
408      zqs = zqs * zcor
409    ELSE
410      IF (t(i, kb)<t_coup) THEN
411        zqs = qsats(t(i, kb)) / pplay(i, kb)
412      ELSE
413        zqs = qsatl(t(i, kb)) / pplay(i, kb)
414      END IF
415    END IF
416    zcll = cloib * zdthmin(i) + cloic
417    zcll = min(1.0, max(0.0, zcll))
418    zrhb = q(i, kb) / zqs
419    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll * (1. - (crhl - zrhb) * cloid)
420    zcll = min(1.0, max(0.0, zcll))
421    diafra(i, kb) = max(diafra(i, kb), zcll)
422    dialiq(i, kb) = diafra(i, kb) * rgammas * zqs
423  END DO
424
425END SUBROUTINE diagcld2
Note: See TracBrowser for help on using the repository browser.