source: LMDZ4/trunk/libf/phytherm/nuage.F @ 1068

Last change on this file since 1068 was 814, checked in by Laurent Fairhead, 17 years ago

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