source: LMDZ6/branches/contrails/libf/phylmd/nuage.f90 @ 5440

Last change on this file since 5440 was 5305, checked in by abarral, 7 weeks ago

Turn YOMCST2.h.h into module
  • 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 5305 2024-10-30 18:29:21Z 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 clesphys_mod_h
7  USE yomcst_mod_h
8  USE dimphy
9  USE lmdz_lscp_tools, only: icefrac_lscp
10  USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
11  USE lmdz_lscp_ini, only : iflag_t_glace
12  USE phys_local_var_mod, ONLY: ptconv
13  USE nuage_params_mod_h
14  IMPLICIT NONE
15  ! ======================================================================
16  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
17  ! Objet: Calculer epaisseur optique et emmissivite des nuages
18  ! ======================================================================
19  ! Arguments:
20  ! t-------input-R-temperature
21  ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
22  ! picefra--inout-R-fraction de glace dans les nuages (-)
23  ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
24  ! ok_aie--input-L-apply aerosol indirect effect or not
25  ! mass_solu_aero-----input-R-total mass concentration for all soluble
26  ! aerosols[ug/m^3]
27  ! mass_solu_aero_pi--input-R-dito, pre-industrial value
28  ! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
29  ! bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
30
31  ! cldtaupi-output-R-pre-industrial value of cloud optical thickness,
32  ! needed for the diagnostics of the aerosol indirect
33  ! radiative forcing (see radlwsw)
34  ! re------output-R-Cloud droplet effective radius multiplied by fl [um]
35  ! fl------output-R-Denominator to re, introduced to avoid problems in
36  ! the averaging of the output. fl is the fraction of liquid
37  ! water clouds within a grid cell
38
39  ! pcltau--output-R-epaisseur optique des nuages
40  ! pclemi--output-R-emissivite des nuages (0 a 1)
41  ! ======================================================================
42
43  REAL paprs(klon, klev+1), pplay(klon, klev)
44  REAL t(klon, klev)
45
46  REAL pclc(klon, klev)
47  REAL pqlwp(klon, klev), picefra(klon,klev)
48  REAL pcltau(klon, klev), pclemi(klon, klev)
49
50  REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
51  REAL distcltop(klon,klev)
52  REAL temp_cltop(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),temp_cltop(:,k),zfice(:),dzfice(:))
116        ELSE
117            CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:))
118
119        ENDIF
120
121    IF (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        IF (.NOT. ptconv(i,k)) THEN
127           zfice(i)=picefra(i,k)
128        ENDIF
129        ENDDO
130    ENDIF
131
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
257  RETURN
258END SUBROUTINE nuage
259SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq)
260  USE dimphy
261  USE yomcst_mod_h
262IMPLICIT 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
272
273  ! Arguments d'entree:
274  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
275  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
276  REAL t(klon, klev) ! temperature (K)
277  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
278  REAL rain(klon) ! pluie convective (kg/m2/s)
279  REAL snow(klon) ! neige convective (kg/m2/s)
280  INTEGER ktop(klon) ! sommet de la convection
281  INTEGER kbot(klon) ! bas de la convection
282
283  ! Arguments de sortie:
284  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
285  REAL dialiq(klon, klev) ! eau liquide nuageuse
286
287  ! Constantes ajustables:
288  REAL canva, canvb, canvh
289  PARAMETER (canva=2.0, canvb=0.3, canvh=0.4)
290  REAL cca, ccb, ccc
291  PARAMETER (cca=0.125, ccb=1.5, ccc=0.8)
292  REAL ccfct, ccscal
293  PARAMETER (ccfct=0.400)
294  PARAMETER (ccscal=1.0E+11)
295  REAL cetahb, cetamb
296  PARAMETER (cetahb=0.45, cetamb=0.80)
297  REAL cclwmr
298  PARAMETER (cclwmr=1.E-04)
299  REAL zepscr
300  PARAMETER (zepscr=1.0E-10)
301
302  ! Variables locales:
303  INTEGER i, k
304  REAL zcc(klon)
305
306  ! Initialisation:
307
308  DO k = 1, klev
309    DO i = 1, klon
310      diafra(i, k) = 0.0
311      dialiq(i, k) = 0.0
312    END DO
313  END DO
314
315  DO i = 1, klon ! Calculer la fraction nuageuse
316    zcc(i) = 0.0
317    IF ((rain(i)+snow(i))>0.) THEN
318      zcc(i) = cca*log(max(zepscr,(rain(i)+snow(i))*ccscal)) - ccb
319      zcc(i) = min(ccc, max(0.0,zcc(i)))
320    END IF
321  END DO
322
323  DO i = 1, klon ! pour traiter les enclumes
324    diafra(i, ktop(i)) = max(diafra(i,ktop(i)), zcc(i)*ccfct)
325    IF ((zcc(i)>=canvh) .AND. (pplay(i,ktop(i))<=cetahb*paprs(i, &
326      1))) diafra(i, ktop(i)) = max(diafra(i,ktop(i)), max(zcc( &
327      i)*ccfct,canva*(zcc(i)-canvb)))
328    dialiq(i, ktop(i)) = cclwmr*diafra(i, ktop(i))
329  END DO
330
331  DO k = 1, klev ! nuages convectifs (sauf enclumes)
332    DO i = 1, klon
333      IF (k<ktop(i) .AND. k>=kbot(i)) THEN
334        diafra(i, k) = max(diafra(i,k), zcc(i)*ccfct)
335        dialiq(i, k) = cclwmr*diafra(i, k)
336      END IF
337    END DO
338  END DO
339
340  RETURN
341END SUBROUTINE diagcld1
342SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
343  USE dimphy
344  USE yomcst_mod_h
345  USE yoethf_mod_h
346IMPLICIT NONE
347
348
349
350  ! Arguments d'entree:
351  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
352  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
353  REAL t(klon, klev) ! temperature (K)
354  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
355
356  ! Arguments de sortie:
357  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
358  REAL dialiq(klon, klev) ! eau liquide nuageuse
359
360  REAL cetamb
361  PARAMETER (cetamb=0.80)
362  REAL cloia, cloib, cloic, cloid
363  PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
364  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
365  REAL rgammas
366  PARAMETER (rgammas=0.05)
367  REAL crhl
368  PARAMETER (crhl=0.15)
369  ! cc      PARAMETER (CRHL=0.70)
370  REAL t_coup
371  PARAMETER (t_coup=234.0)
372
373  ! Variables locales:
374  INTEGER i, k, kb, invb(klon)
375  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
376  REAL zdelta, zcor
377
378  ! Fonctions thermodynamiques:
379  include "FCTTRE.h"
380
381  ! Initialisation:
382
383  DO k = 1, klev
384    DO i = 1, klon
385      diafra(i, k) = 0.0
386      dialiq(i, k) = 0.0
387    END DO
388  END DO
389
390  DO i = 1, klon
391    invb(i) = klev
392    zdthmin(i) = 0.0
393  END DO
394
395  DO k = 2, klev/2 - 1
396    DO i = 1, klon
397      zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &
398        rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)
399      zdthdp = zdthdp*cloia
400      IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN
401        zdthmin(i) = zdthdp
402        invb(i) = k
403      END IF
404    END DO
405  END DO
406
407  DO i = 1, klon
408    kb = invb(i)
409    IF (thermcep) THEN
410      zdelta = max(0., sign(1.,rtt-t(i,kb)))
411      zqs = r2es*foeew(t(i,kb), zdelta)/pplay(i, kb)
412      zqs = min(0.5, zqs)
413      zcor = 1./(1.-retv*zqs)
414      zqs = zqs*zcor
415    ELSE
416      IF (t(i,kb)<t_coup) THEN
417        zqs = qsats(t(i,kb))/pplay(i, kb)
418      ELSE
419        zqs = qsatl(t(i,kb))/pplay(i, kb)
420      END IF
421    END IF
422    zcll = cloib*zdthmin(i) + cloic
423    zcll = min(1.0, max(0.0,zcll))
424    zrhb = q(i, kb)/zqs
425    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid)
426    zcll = min(1.0, max(0.0,zcll))
427    diafra(i, kb) = max(diafra(i,kb), zcll)
428    dialiq(i, kb) = diafra(i, kb)*rgammas*zqs
429  END DO
430
431  RETURN
432END SUBROUTINE diagcld2
Note: See TracBrowser for help on using the repository browser.