Context Navigation

nuage.F @ 1267

Last change on this file since 1267 was 1183, checked in by jghattas, 15 years ago

Ajoute test de coherence entre les longitudes et latitudes des fichiers de forcage des aerosols et le model.

Corection pour la variable mass_solu_aero : changement du nom, commentaire et contenu(il y avait l'aerosol dust en trop)

Changement du nom de module car dans INCA le meme nom existe deja : aerosol_mod => aero_mod

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 13.7 KB

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

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

Download in other formats:

Original Format