source: LMDZ5/trunk/libf/phylmd/nuage.F90 @ 2077

Last change on this file since 2077 was 2077, checked in by fhourdin, 10 years ago

Modification relative à la phase mixte des nuages :

  1. on fait tendre t_glace_min vers t_glace_max (en principe 0°C) linéairerment entre p/ps = 0.8 et p/ps=1.
  2. Passage de tableau à la routine icefrac_lsc en remplacement d'une fonction scalaire.
  3. Changement de nom pour le module (le module icefrac_lsc_mod.F90 contenant maintenant la routine icefrac_lsc).

Mofication concerning the mixte phase of clouds

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