source: LMDZ5/branches/testing/libf/phylmd/nuage.F90 @ 2176

Last change on this file since 2176 was 2160, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
  • 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 2160 2014-11-28 15:36:29Z jescribano $
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 ! computes ice fraction (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 icefrac_lsc_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),zfice(:))
111     ENDIF
112
113    DO i = 1, klon
114      rad_chaud = rad_chau1
115      IF (k<=3) rad_chaud = rad_chau2
116
117      pclc(i, k) = max(pclc(i,k), seuil_neb)
118      zflwp = 1000.*pqlwp(i, k)/rg/pclc(i, k)*(paprs(i,k)-paprs(i,k+1))
119
120      IF (ok_aie) THEN
121          ! Formula "D" of Boucher and Lohmann, Tellus, 1995
122          !
123        cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
124          1.E-4))/log(10.))*1.E6 !-m-3
125          ! Cloud droplet number concentration (CDNC) is restricted
126          ! to be within [20, 1000 cm^3]
127          !
128        cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k)))
129        cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
130          1.E-4))/log(10.))*1.E6 !-m-3
131        cdnc_pi(i, k) = min(1000.E6, max(20.E6,cdnc_pi(i,k)))
132          !
133          !
134          ! air density: pplay(i,k) / (RD * zT(i,k))
135          ! factor 1.1: derive effective radius from volume-mean radius
136          ! factor 1000 is the water density
137          ! _chaud means that this is the CDR for liquid water clouds
138          !
139        rad_chaud = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi*1000. &
140          *cdnc(i,k)))**(1./3.)
141          !
142          ! Convert to um. CDR shall be at least 3 um.
143          !
144        rad_chaud = max(rad_chaud*1.E6, 3.)
145
146          ! For output diagnostics
147          !
148          ! Cloud droplet effective radius [um]
149          !
150          ! we multiply here with f * xl (fraction of liquid water
151          ! clouds in the grid cell) to avoid problems in the
152          ! averaging of the output.
153          ! In the output of IOIPSL, derive the real cloud droplet
154          ! effective radius as re/fl
155          !
156        fl(i, k) = pclc(i, k)*(1.-zfice(i))
157        re(i, k) = rad_chaud*fl(i, k)
158
159          ! Pre-industrial cloud opt thickness
160          !
161          ! "radius" is calculated as rad_chaud above (plus the
162          ! ice cloud contribution) but using cdnc_pi instead of
163          ! cdnc.
164        radius = max(1.1E6*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi* &
165          1000.*cdnc_pi(i,k)))**(1./3.), 3.)*(1.-zfice(i)) + rad_froid*zfice(i)
166        cldtaupi(i, k) = 3.0/2.0*zflwp/radius
167      END IF ! ok_aie
168
169      radius = rad_chaud*(1.-zfice(i)) + rad_froid*zfice(i)
170      coef = coef_chau*(1.-zfice(i)) + coef_froi*zfice(i)
171      pcltau(i, k) = 3.0/2.0*zflwp/radius
172      pclemi(i, k) = 1.0 - exp(-coef*zflwp)
173      lo = (pclc(i,k)<=seuil_neb)
174      IF (lo) pclc(i, k) = 0.0
175      IF (lo) pcltau(i, k) = 0.0
176      IF (lo) pclemi(i, k) = 0.0
177
178      IF (.NOT. ok_aie) cldtaupi(i, k) = pcltau(i, k)
179    END DO
180  END DO
181  ! cc      DO k = 1, klev
182  ! cc      DO i = 1, klon
183  ! cc         t(i,k) = t(i,k)
184  ! cc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
185  ! cc         lo = pclc(i,k) .GT. (2.*1.e-5)
186  ! cc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
187  ! cc     .          /(rg*pclc(i,k))
188  ! cc         zradef = 10.0 + (1.-sigs(k))*45.0
189  ! cc         pcltau(i,k) = 1.5 * zflwp / zradef
190  ! cc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
191  ! cc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
192  ! cc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
193  ! cc         if (.NOT.lo) pclc(i,k) = 0.0
194  ! cc         if (.NOT.lo) pcltau(i,k) = 0.0
195  ! cc         if (.NOT.lo) pclemi(i,k) = 0.0
196  ! cc      ENDDO
197  ! cc      ENDDO
198  ! ccccc      print*, 'pas de nuage dans le rayonnement'
199  ! ccccc      DO k = 1, klev
200  ! ccccc      DO i = 1, klon
201  ! ccccc         pclc(i,k) = 0.0
202  ! ccccc         pcltau(i,k) = 0.0
203  ! ccccc         pclemi(i,k) = 0.0
204  ! ccccc      ENDDO
205  ! ccccc      ENDDO
206
207  ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
208
209  DO i = 1, klon
210    pct(i) = 1.0
211    pch(i) = 1.0
212    pcm(i) = 1.0
213    pcl(i) = 1.0
214    pctlwp(i) = 0.0
215  END DO
216
217  DO k = klev, 1, -1
218    DO i = 1, klon
219      pctlwp(i) = pctlwp(i) + pqlwp(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
220      pct(i) = pct(i)*(1.0-pclc(i,k))
221      IF (pplay(i,k)<=cetahb*paprs(i,1)) pch(i) = pch(i)*(1.0-pclc(i,k))
222      IF (pplay(i,k)>cetahb*paprs(i,1) .AND. pplay(i,k)<=cetamb*paprs(i,1)) &
223        pcm(i) = pcm(i)*(1.0-pclc(i,k))
224      IF (pplay(i,k)>cetamb*paprs(i,1)) pcl(i) = pcl(i)*(1.0-pclc(i,k))
225    END DO
226  END DO
227
228  DO i = 1, klon
229    pct(i) = 1. - pct(i)
230    pch(i) = 1. - pch(i)
231    pcm(i) = 1. - pcm(i)
232    pcl(i) = 1. - pcl(i)
233  END DO
234
235  RETURN
236END SUBROUTINE nuage
237SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq)
238  USE dimphy
239  IMPLICIT NONE
240
241  ! Laurent Li (LMD/CNRS), le 12 octobre 1998
242  ! (adaptation du code ECMWF)
243
244  ! Dans certains cas, le schema pronostique des nuages n'est
245  ! pas suffisament performant. On a donc besoin de diagnostiquer
246  ! ces nuages. Je dois avouer que c'est une frustration.
247
248  ! ym#include "dimensions.h"
249  ! ym#include "dimphy.h"
250  include "YOMCST.h"
251
252  ! Arguments d'entree:
253  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
254  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
255  REAL t(klon, klev) ! temperature (K)
256  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
257  REAL rain(klon) ! pluie convective (kg/m2/s)
258  REAL snow(klon) ! neige convective (kg/m2/s)
259  INTEGER ktop(klon) ! sommet de la convection
260  INTEGER kbot(klon) ! bas de la convection
261
262  ! Arguments de sortie:
263  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
264  REAL dialiq(klon, klev) ! eau liquide nuageuse
265
266  ! Constantes ajustables:
267  REAL canva, canvb, canvh
268  PARAMETER (canva=2.0, canvb=0.3, canvh=0.4)
269  REAL cca, ccb, ccc
270  PARAMETER (cca=0.125, ccb=1.5, ccc=0.8)
271  REAL ccfct, ccscal
272  PARAMETER (ccfct=0.400)
273  PARAMETER (ccscal=1.0E+11)
274  REAL cetahb, cetamb
275  PARAMETER (cetahb=0.45, cetamb=0.80)
276  REAL cclwmr
277  PARAMETER (cclwmr=1.E-04)
278  REAL zepscr
279  PARAMETER (zepscr=1.0E-10)
280
281  ! Variables locales:
282  INTEGER i, k
283  REAL zcc(klon)
284
285  ! Initialisation:
286
287  DO k = 1, klev
288    DO i = 1, klon
289      diafra(i, k) = 0.0
290      dialiq(i, k) = 0.0
291    END DO
292  END DO
293
294  DO i = 1, klon ! Calculer la fraction nuageuse
295    zcc(i) = 0.0
296    IF ((rain(i)+snow(i))>0.) THEN
297      zcc(i) = cca*log(max(zepscr,(rain(i)+snow(i))*ccscal)) - ccb
298      zcc(i) = min(ccc, max(0.0,zcc(i)))
299    END IF
300  END DO
301
302  DO i = 1, klon ! pour traiter les enclumes
303    diafra(i, ktop(i)) = max(diafra(i,ktop(i)), zcc(i)*ccfct)
304    IF ((zcc(i)>=canvh) .AND. (pplay(i,ktop(i))<=cetahb*paprs(i, &
305      1))) diafra(i, ktop(i)) = max(diafra(i,ktop(i)), max(zcc( &
306      i)*ccfct,canva*(zcc(i)-canvb)))
307    dialiq(i, ktop(i)) = cclwmr*diafra(i, ktop(i))
308  END DO
309
310  DO k = 1, klev ! nuages convectifs (sauf enclumes)
311    DO i = 1, klon
312      IF (k<ktop(i) .AND. k>=kbot(i)) THEN
313        diafra(i, k) = max(diafra(i,k), zcc(i)*ccfct)
314        dialiq(i, k) = cclwmr*diafra(i, k)
315      END IF
316    END DO
317  END DO
318
319  RETURN
320END SUBROUTINE diagcld1
321SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
322  USE dimphy
323  IMPLICIT NONE
324
325  ! ym#include "dimensions.h"
326  ! ym#include "dimphy.h"
327  include "YOMCST.h"
328
329  ! Arguments d'entree:
330  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
331  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
332  REAL t(klon, klev) ! temperature (K)
333  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
334
335  ! Arguments de sortie:
336  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
337  REAL dialiq(klon, klev) ! eau liquide nuageuse
338
339  REAL cetamb
340  PARAMETER (cetamb=0.80)
341  REAL cloia, cloib, cloic, cloid
342  PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
343  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
344  REAL rgammas
345  PARAMETER (rgammas=0.05)
346  REAL crhl
347  PARAMETER (crhl=0.15)
348  ! cc      PARAMETER (CRHL=0.70)
349  REAL t_coup
350  PARAMETER (t_coup=234.0)
351
352  ! Variables locales:
353  INTEGER i, k, kb, invb(klon)
354  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
355  REAL zdelta, zcor
356
357  ! Fonctions thermodynamiques:
358  include "YOETHF.h"
359  include "FCTTRE.h"
360
361  ! Initialisation:
362
363  DO k = 1, klev
364    DO i = 1, klon
365      diafra(i, k) = 0.0
366      dialiq(i, k) = 0.0
367    END DO
368  END DO
369
370  DO i = 1, klon
371    invb(i) = klev
372    zdthmin(i) = 0.0
373  END DO
374
375  DO k = 2, klev/2 - 1
376    DO i = 1, klon
377      zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &
378        rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)
379      zdthdp = zdthdp*cloia
380      IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN
381        zdthmin(i) = zdthdp
382        invb(i) = k
383      END IF
384    END DO
385  END DO
386
387  DO i = 1, klon
388    kb = invb(i)
389    IF (thermcep) THEN
390      zdelta = max(0., sign(1.,rtt-t(i,kb)))
391      zqs = r2es*foeew(t(i,kb), zdelta)/pplay(i, kb)
392      zqs = min(0.5, zqs)
393      zcor = 1./(1.-retv*zqs)
394      zqs = zqs*zcor
395    ELSE
396      IF (t(i,kb)<t_coup) THEN
397        zqs = qsats(t(i,kb))/pplay(i, kb)
398      ELSE
399        zqs = qsatl(t(i,kb))/pplay(i, kb)
400      END IF
401    END IF
402    zcll = cloib*zdthmin(i) + cloic
403    zcll = min(1.0, max(0.0,zcll))
404    zrhb = q(i, kb)/zqs
405    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid)
406    zcll = min(1.0, max(0.0,zcll))
407    diafra(i, kb) = max(diafra(i,kb), zcll)
408    dialiq(i, kb) = diafra(i, kb)*rgammas*zqs
409  END DO
410
411  RETURN
412END SUBROUTINE diagcld2
Note: See TracBrowser for help on using the repository browser.