source: LMDZ6/trunk/libf/phylmd/nuage.F90 @ 4639

Last change on this file since 4639 was 4639, checked in by evignon, 15 months ago

modification de Lea pour inclure une dependance de la phase a la temperature du sommet du nuage
  • 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.7 KB
Line 
1! $Id: nuage.F90 4639 2023-07-24 14:46:18Z evignon $
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 lscp_tools_mod, only: icefrac_lscp
8  USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
9  USE lscp_ini_mod, only : iflag_t_glace
10  USE phys_local_var_mod, ONLY: ptconv
11  IMPLICIT NONE
12  ! ======================================================================
13  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
14  ! Objet: Calculer epaisseur optique et emmissivite des nuages
15  ! ======================================================================
16  ! Arguments:
17  ! t-------input-R-temperature
18  ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
19  ! picefra--inout-R-fraction de glace dans les nuages (-)
20  ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
21  ! ok_aie--input-L-apply aerosol indirect effect or not
22  ! mass_solu_aero-----input-R-total mass concentration for all soluble
23  ! aerosols[ug/m^3]
24  ! mass_solu_aero_pi--input-R-dito, pre-industrial value
25  ! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
26  ! bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
27
28  ! cldtaupi-output-R-pre-industrial value of cloud optical thickness,
29  ! needed for the diagnostics of the aerosol indirect
30  ! radiative forcing (see radlwsw)
31  ! re------output-R-Cloud droplet effective radius multiplied by fl [um]
32  ! fl------output-R-Denominator to re, introduced to avoid problems in
33  ! the averaging of the output. fl is the fraction of liquid
34  ! water clouds within a grid cell
35
36  ! pcltau--output-R-epaisseur optique des nuages
37  ! pclemi--output-R-emissivite des nuages (0 a 1)
38  ! ======================================================================
39
40  include "YOMCST.h"
41  include "nuage.h" ! JBM 3/14
42  include "clesphys.h"
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.EQ.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 ((.NOT. ptconv(i,k)) .AND. 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            zfice(i)=picefra(i,k)
128            ENDDO
129        ENDIF
130
131
132     ENDIF
133
134    DO i = 1, klon
135      rad_chaud = rad_chau1
136      IF (k<=3) rad_chaud = rad_chau2
137
138      pclc(i, k) = max(pclc(i,k), seuil_neb)
139      zflwp = 1000.*pqlwp(i, k)/rg/pclc(i, k)*(paprs(i,k)-paprs(i,k+1))
140
141      IF (ok_aie) THEN
142          ! Formula "D" of Boucher and Lohmann, Tellus, 1995
143          !
144        cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
145          1.E-4))/log(10.))*1.E6 !-m-3
146          ! Cloud droplet number concentration (CDNC) is restricted
147          ! to be within [20, 1000 cm^3]
148          !
149        cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k)))
150        cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
151          1.E-4))/log(10.))*1.E6 !-m-3
152        cdnc_pi(i, k) = min(1000.E6, max(20.E6,cdnc_pi(i,k)))
153          !
154          !
155          ! air density: pplay(i,k) / (RD * zT(i,k))
156          ! factor 1.1: derive effective radius from volume-mean radius
157          ! factor 1000 is the water density
158          ! _chaud means that this is the CDR for liquid water clouds
159          !
160        rad_chaud = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi*1000. &
161          *cdnc(i,k)))**(1./3.)
162          !
163          ! Convert to um. CDR shall be at least 3 um.
164          !
165        rad_chaud = max(rad_chaud*1.E6, 3.)
166
167          ! For output diagnostics
168          !
169          ! Cloud droplet effective radius [um]
170          !
171          ! we multiply here with f * xl (fraction of liquid water
172          ! clouds in the grid cell) to avoid problems in the
173          ! averaging of the output.
174          ! In the output of IOIPSL, derive the real cloud droplet
175          ! effective radius as re/fl
176          !
177        fl(i, k) = pclc(i, k)*(1.-zfice(i))
178        re(i, k) = rad_chaud*fl(i, k)
179
180          ! Pre-industrial cloud opt thickness
181          !
182          ! "radius" is calculated as rad_chaud above (plus the
183          ! ice cloud contribution) but using cdnc_pi instead of
184          ! cdnc.
185        radius = max(1.1E6*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi* &
186          1000.*cdnc_pi(i,k)))**(1./3.), 3.)*(1.-zfice(i)) + rad_froid*zfice(i)
187        cldtaupi(i, k) = 3.0/2.0*zflwp/radius
188      END IF ! ok_aie
189
190      radius = rad_chaud*(1.-zfice(i)) + rad_froid*zfice(i)
191      coef = coef_chau*(1.-zfice(i)) + coef_froi*zfice(i)
192      pcltau(i, k) = 3.0/2.0*zflwp/radius
193      pclemi(i, k) = 1.0 - exp(-coef*zflwp)
194      lo = (pclc(i,k)<=seuil_neb)
195      IF (lo) pclc(i, k) = 0.0
196      IF (lo) pcltau(i, k) = 0.0
197      IF (lo) pclemi(i, k) = 0.0
198
199      IF (.NOT. ok_aie) cldtaupi(i, k) = pcltau(i, k)
200    END DO
201  END DO
202  ! cc      DO k = 1, klev
203  ! cc      DO i = 1, klon
204  ! cc         t(i,k) = t(i,k)
205  ! cc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
206  ! cc         lo = pclc(i,k) .GT. (2.*1.e-5)
207  ! cc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
208  ! cc     .          /(rg*pclc(i,k))
209  ! cc         zradef = 10.0 + (1.-sigs(k))*45.0
210  ! cc         pcltau(i,k) = 1.5 * zflwp / zradef
211  ! cc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
212  ! cc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
213  ! cc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
214  ! cc         if (.NOT.lo) pclc(i,k) = 0.0
215  ! cc         if (.NOT.lo) pcltau(i,k) = 0.0
216  ! cc         if (.NOT.lo) pclemi(i,k) = 0.0
217  ! cc      ENDDO
218  ! cc      ENDDO
219  ! ccccc      print*, 'pas de nuage dans le rayonnement'
220  ! ccccc      DO k = 1, klev
221  ! ccccc      DO i = 1, klon
222  ! ccccc         pclc(i,k) = 0.0
223  ! ccccc         pcltau(i,k) = 0.0
224  ! ccccc         pclemi(i,k) = 0.0
225  ! ccccc      ENDDO
226  ! ccccc      ENDDO
227
228  ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
229
230  DO i = 1, klon
231    pct(i) = 1.0
232    pch(i) = 1.0
233    pcm(i) = 1.0
234    pcl(i) = 1.0
235    pctlwp(i) = 0.0
236  END DO
237
238  DO k = klev, 1, -1
239    DO i = 1, klon
240      pctlwp(i) = pctlwp(i) + pqlwp(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
241      pct(i) = pct(i)*(1.0-pclc(i,k))
242      IF (pplay(i,k)<=cetahb*paprs(i,1)) pch(i) = pch(i)*(1.0-pclc(i,k))
243      IF (pplay(i,k)>cetahb*paprs(i,1) .AND. pplay(i,k)<=cetamb*paprs(i,1)) &
244        pcm(i) = pcm(i)*(1.0-pclc(i,k))
245      IF (pplay(i,k)>cetamb*paprs(i,1)) pcl(i) = pcl(i)*(1.0-pclc(i,k))
246    END DO
247  END DO
248
249  DO i = 1, klon
250    pct(i) = 1. - pct(i)
251    pch(i) = 1. - pch(i)
252    pcm(i) = 1. - pcm(i)
253    pcl(i) = 1. - pcl(i)
254  END DO
255
256  RETURN
257END SUBROUTINE nuage
258SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq)
259  USE dimphy
260  IMPLICIT NONE
261
262  ! Laurent Li (LMD/CNRS), le 12 octobre 1998
263  ! (adaptation du code ECMWF)
264
265  ! Dans certains cas, le schema pronostique des nuages n'est
266  ! pas suffisament performant. On a donc besoin de diagnostiquer
267  ! ces nuages. Je dois avouer que c'est une frustration.
268
269  include "YOMCST.h"
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
338  RETURN
339END SUBROUTINE diagcld1
340SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
341  USE dimphy
342  IMPLICIT NONE
343
344  include "YOMCST.h"
345
346  ! Arguments d'entree:
347  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
348  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
349  REAL t(klon, klev) ! temperature (K)
350  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
351
352  ! Arguments de sortie:
353  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
354  REAL dialiq(klon, klev) ! eau liquide nuageuse
355
356  REAL cetamb
357  PARAMETER (cetamb=0.80)
358  REAL cloia, cloib, cloic, cloid
359  PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
360  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
361  REAL rgammas
362  PARAMETER (rgammas=0.05)
363  REAL crhl
364  PARAMETER (crhl=0.15)
365  ! cc      PARAMETER (CRHL=0.70)
366  REAL t_coup
367  PARAMETER (t_coup=234.0)
368
369  ! Variables locales:
370  INTEGER i, k, kb, invb(klon)
371  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
372  REAL zdelta, zcor
373
374  ! Fonctions thermodynamiques:
375  include "YOETHF.h"
376  include "FCTTRE.h"
377
378  ! Initialisation:
379
380  DO k = 1, klev
381    DO i = 1, klon
382      diafra(i, k) = 0.0
383      dialiq(i, k) = 0.0
384    END DO
385  END DO
386
387  DO i = 1, klon
388    invb(i) = klev
389    zdthmin(i) = 0.0
390  END DO
391
392  DO k = 2, klev/2 - 1
393    DO i = 1, klon
394      zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &
395        rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)
396      zdthdp = zdthdp*cloia
397      IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN
398        zdthmin(i) = zdthdp
399        invb(i) = k
400      END IF
401    END DO
402  END DO
403
404  DO i = 1, klon
405    kb = invb(i)
406    IF (thermcep) THEN
407      zdelta = max(0., sign(1.,rtt-t(i,kb)))
408      zqs = r2es*foeew(t(i,kb), zdelta)/pplay(i, kb)
409      zqs = min(0.5, zqs)
410      zcor = 1./(1.-retv*zqs)
411      zqs = zqs*zcor
412    ELSE
413      IF (t(i,kb)<t_coup) THEN
414        zqs = qsats(t(i,kb))/pplay(i, kb)
415      ELSE
416        zqs = qsatl(t(i,kb))/pplay(i, kb)
417      END IF
418    END IF
419    zcll = cloib*zdthmin(i) + cloic
420    zcll = min(1.0, max(0.0,zcll))
421    zrhb = q(i, kb)/zqs
422    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid)
423    zcll = min(1.0, max(0.0,zcll))
424    diafra(i, kb) = max(diafra(i,kb), zcll)
425    dialiq(i, kb) = diafra(i, kb)*rgammas*zqs
426  END DO
427
428  RETURN
429END SUBROUTINE diagcld2
Note: See TracBrowser for help on using the repository browser.