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

Last change on this file since 4283 was 3999, checked in by evignon, 3 years ago

commission de la nouvelle routine de condensation
grande echelle simplifiee (lscp, version epuree de fistilp)
et du schema de nuages de phase mixte (en developpement)

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