source: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/nuage.F90 @ 5408

Last change on this file since 5408 was 4669, checked in by Laurent Fairhead, 16 months ago

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