Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

nuage.F

Last change on this file 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

Line
1	SUBROUTINE nuage (paprs, pplay,
2	. t, pqlwp, pclc, pcltau, pclemi,
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)
8	IMPLICIT none
9	c======================================================================
10	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
11	c Objet: Calculer epaisseur optique et emmissivite des nuages
12	c======================================================================
13	c Arguments:
14	c t-------input-R-temperature
15	c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
16	c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
17	c ok_aie--input-L-apply aerosol indirect effect or not
18	c sulfate-input-R-sulfate aerosol mass concentration [um/m^3]
19	c sulfate_pi-input-R-dito, pre-industrial value
20	c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
21	c bl95_b1-input-R-a parameter, may be varied for tests ( -"- )
22	c
23	c cldtaupi-output-R-pre-industrial value of cloud optical thickness,
24	c needed for the diagnostics of the aerosol indirect
25	c radiative forcing (see radlwsw)
26	c re------output-R-Cloud droplet effective radius multiplied by fl [um]
27	c fl------output-R-Denominator to re, introduced to avoid problems in
28	c the averaging of the output. fl is the fraction of liquid
29	c water clouds within a grid cell
30	c
31	c pcltau--output-R-epaisseur optique des nuages
32	c pclemi--output-R-emissivite des nuages (0 a 1)
33	c======================================================================
34	C
35	#include "YOMCST.h"
36	c
37	#include "dimensions.h"
38	#include "dimphy.h"
39	REAL paprs(klon,klev+1), pplay(klon,klev)
40	REAL t(klon,klev)
41	c
42	REAL pclc(klon,klev)
43	REAL pqlwp(klon,klev)
44	REAL pcltau(klon,klev), pclemi(klon,klev)
45	c
46	REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
47	c
48	LOGICAL lo
49	c
50	REAL cetahb, cetamb
51	PARAMETER (cetahb = 0.45, cetamb = 0.80)
52	C
53	INTEGER i, k
54	REAL zflwp, zradef, zfice, zmsac
55	c
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)
58	ccc PARAMETER (rad_chaud=15.0, rad_froid=35.0)
59	c sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0)
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)
66
67	cjq for the aerosol indirect effect
68	cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
69	cjq
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
84	cjq-end
85
86	ccc PARAMETER (nexpo=1)
87	c
88	c Calculer l'epaisseur optique et l'emmissivite des nuages
89	c
90	DO k = 1, klev
91	DO i = 1, klon
92	rad_chaud = rad_chau1
93	IF (k.LE.3) rad_chaud = rad_chau2
94
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
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)/(RDT(i,k)))
148	. / (4./3.RPI1000.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
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
162
163	IF (.NOT.ok_aie) cldtaupi(i,k)=pcltau(i,k)
164	ENDDO
165	ENDDO
166	ccc DO k = 1, klev
167	ccc DO i = 1, klon
168	ccc t(i,k) = t(i,k)
169	ccc pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
170	ccc lo = pclc(i,k) .GT. (2.*1.e-5)
171	ccc zflwp = pqlwp(i,k)1000.(paprs(i,k)-paprs(i,k+1))
172	ccc . /(rg*pclc(i,k))
173	ccc zradef = 10.0 + (1.-sigs(k))*45.0
174	ccc pcltau(i,k) = 1.5 * zflwp / zradef
175	ccc zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
176	ccc zmsac = 0.13(1.0-zfice) + 0.08zfice
177	ccc pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
178	ccc if (.NOT.lo) pclc(i,k) = 0.0
179	ccc if (.NOT.lo) pcltau(i,k) = 0.0
180	ccc if (.NOT.lo) pclemi(i,k) = 0.0
181	ccc ENDDO
182	ccc ENDDO
183	cccccc print*, 'pas de nuage dans le rayonnement'
184	cccccc DO k = 1, klev
185	cccccc DO i = 1, klon
186	cccccc pclc(i,k) = 0.0
187	cccccc pcltau(i,k) = 0.0
188	cccccc pclemi(i,k) = 0.0
189	cccccc ENDDO
190	cccccc ENDDO
191	C
192	C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
193	C
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
201	C
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
216	C
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
223	C
224	RETURN
225	END
226	SUBROUTINE diagcld1(paprs,pplay,rain,snow,kbot,ktop,
227	. diafra,dialiq)
228	IMPLICIT none
229	c
230	c Laurent Li (LMD/CNRS), le 12 octobre 1998
231	c (adaptation du code ECMWF)
232	c
233	c Dans certains cas, le schema pronostique des nuages n'est
234	c pas suffisament performant. On a donc besoin de diagnostiquer
235	c ces nuages. Je dois avouer que c'est une frustration.
236	c
237	#include "dimensions.h"
238	#include "dimphy.h"
239	#include "YOMCST.h"
240	c
241	c 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
250	c
251	c Arguments de sortie:
252	REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
253	REAL dialiq(klon,klev) ! eau liquide nuageuse
254	c
255	c 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)
269	c
270	c Variables locales:
271	INTEGER i, k
272	REAL zcc(klon)
273	c
274	c Initialisation:
275	c
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
282	c
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
290	c
291	DO i = 1, klon ! pour traiter les enclumes
292	diafra(i,ktop(i)) = MAX(diafra(i,ktop(i)),zcc(i)*CCFCT)
293	IF ((zcc(i).GE.CANVH) .AND.
294	. (pplay(i,ktop(i)).LE.CETAHB*paprs(i,1)))
295	. diafra(i,ktop(i)) = MAX(diafra(i,ktop(i)),
296	. MAX(zcc(i)CCFCT,CANVA(zcc(i)-CANVB)))
297	dialiq(i,ktop(i))=CCLWMR*diafra(i,ktop(i))
298	ENDDO
299	c
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
303	diafra(i,k)=MAX(diafra(i,k),zcc(i)*CCFCT)
304	dialiq(i,k)=CCLWMR*diafra(i,k)
305	ENDIF
306	ENDDO
307	ENDDO
308	c
309	RETURN
310	END
311	SUBROUTINE diagcld2(paprs,pplay,t,q, diafra,dialiq)
312	IMPLICIT none
313	c
314	#include "dimensions.h"
315	#include "dimphy.h"
316	#include "YOMCST.h"
317	c
318	c 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)
323	c
324	c Arguments de sortie:
325	REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
326	REAL dialiq(klon,klev) ! eau liquide nuageuse
327	c
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)
332	ccc 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)
337	ccc PARAMETER (CRHL=0.70)
338	REAL t_coup
339	PARAMETER (t_coup=234.0)
340	c
341	c Variables locales:
342	INTEGER i, k, kb, invb(klon)
343	REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
344	REAL zdelta, zcor
345	c
346	c Fonctions thermodynamiques:
347	#include "YOETHF.h"
348	#include "FCTTRE.h"
349	c
350	c Initialisation:
351	c
352	DO k = 1, klev
353	DO i = 1, klon
354	diafra(i,k) = 0.0
355	dialiq(i,k) = 0.0
356	ENDDO
357	ENDDO
358	c
359	DO i = 1, klon
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))
398	diafra(i,kb) = MAX(diafra(i,kb),zcll)
399	dialiq(i,kb)= diafra(i,kb) * RGAMMAS*zqs
400	ENDDO
401	c
402	RETURN
403	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format