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

Last change on this file since 5218 was 5153, checked in by abarral, 5 months ago

Revert FCTTRE to INCLUDE to assess impact of inlining
  • 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: 13.9 KB
Line 
1! $Id: nuage.F90 5153 2024-07-31 16:20:03Z 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
343  USE lmdz_yomcst
344
345  IMPLICIT NONE
346 INCLUDE "FCTTRE.h"
347
348  ! Arguments d'entree:
349  REAL paprs(klon, klev + 1) ! pression (Pa) a inter-couche
350  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
351  REAL t(klon, klev) ! temperature (K)
352  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
353
354  ! Arguments de sortie:
355  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
356  REAL dialiq(klon, klev) ! eau liquide nuageuse
357
358  REAL cetamb
359  PARAMETER (cetamb = 0.80)
360  REAL cloia, cloib, cloic, cloid
361  PARAMETER (cloia = 1.0E+02, cloib = -10.00, cloic = -0.6, cloid = 5.0)
362  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
363  REAL rgammas
364  PARAMETER (rgammas = 0.05)
365  REAL crhl
366  PARAMETER (crhl = 0.15)
367  ! cc      PARAMETER (CRHL=0.70)
368  REAL t_coup
369  PARAMETER (t_coup = 234.0)
370
371  ! Variables locales:
372  INTEGER i, k, kb, invb(klon)
373  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
374  REAL zdelta, zcor
375
376  ! Initialisation:
377
378  DO k = 1, klev
379    DO i = 1, klon
380      diafra(i, k) = 0.0
381      dialiq(i, k) = 0.0
382    END DO
383  END DO
384
385  DO i = 1, klon
386    invb(i) = klev
387    zdthmin(i) = 0.0
388  END DO
389
390  DO k = 2, klev / 2 - 1
391    DO i = 1, klon
392      zdthdp = (t(i, k) - t(i, k + 1)) / (pplay(i, k) - pplay(i, k + 1)) - &
393              rd * 0.5 * (t(i, k) + t(i, k + 1)) / rcpd / paprs(i, k + 1)
394      zdthdp = zdthdp * cloia
395      IF (pplay(i, k)>cetamb * paprs(i, 1) .AND. zdthdp<zdthmin(i)) THEN
396        zdthmin(i) = zdthdp
397        invb(i) = k
398      END IF
399    END DO
400  END DO
401
402  DO i = 1, klon
403    kb = invb(i)
404    IF (thermcep) THEN
405      zdelta = max(0., sign(1., rtt - t(i, kb)))
406      zqs = r2es * foeew(t(i, kb), zdelta) / pplay(i, kb)
407      zqs = min(0.5, zqs)
408      zcor = 1. / (1. - retv * zqs)
409      zqs = zqs * zcor
410    ELSE
411      IF (t(i, kb)<t_coup) THEN
412        zqs = qsats(t(i, kb)) / pplay(i, kb)
413      ELSE
414        zqs = qsatl(t(i, kb)) / pplay(i, kb)
415      END IF
416    END IF
417    zcll = cloib * zdthmin(i) + cloic
418    zcll = min(1.0, max(0.0, zcll))
419    zrhb = q(i, kb) / zqs
420    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll * (1. - (crhl - zrhb) * cloid)
421    zcll = min(1.0, max(0.0, zcll))
422    diafra(i, kb) = max(diafra(i, kb), zcll)
423    dialiq(i, kb) = diafra(i, kb) * rgammas * zqs
424  END DO
425
426END SUBROUTINE diagcld2
Note: See TracBrowser for help on using the repository browser.