source: LMDZ.3.3/branches/rel-LF/libf/phylmd/nuage.F @ 2168

Last change on this file since 2168 was 517, checked in by lmdzadmin, 21 years ago

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